Structured FISTA for Image Restoration

TL;DR

This paper introduces sFISTA, an efficient algorithm for large-scale image restoration that exploits matrix structures and Kronecker product approximations to accelerate computation while maintaining accuracy.

Contribution

The paper develops sFISTA, a novel structured FISTA variant that leverages matrix structures and Kronecker approximations for faster large-scale image restoration.

Findings

sFISTA achieves comparable accuracy to FISTA.

sFISTA significantly reduces computational complexity.

Numerical results confirm efficiency improvements.

Abstract

In this paper, we propose an efficient numerical scheme for solving some large scale ill-posed linear inverse problems arising from image restoration. In order to accelerate the computation, two different hidden structures are exploited. First, the coefficient matrix is approximated as the sum of a small number of Kronecker products. This procedure not only introduces one more level of parallelism into the computation but also enables the usage of computationally intensive matrix-matrix multiplications in the subsequent optimization procedure. We then derive the corresponding Tikhonov regularized minimization model and extend the fast iterative shrinkage-thresholding algorithm (FISTA) to solve the resulting optimization problem. Since the matrices appearing in the Kronecker product approximation are all structured matrices (Toeplitz, Hankel, etc.), we can further exploit their fast matrix-vector multiplication algorithms at each iteration. The proposed algorithm is thus called structured fast iterative shrinkage-thresholding algorithm (sFISTA). In particular, we show that the approximation error introduced by sFISTA is well under control and sFISTA can reach the same image restoration accuracy level as FISTA. Finally, both the theoretical complexity analysis and some numerical results are provided to demonstrate the efficiency of sFISTA.

Tables

Table 1. Table 1 : Numerical results for FISTA and sFISTA for ‘hst’

		FISTA		sFISTA		sFISTA		sFISTA		sFISTA		sFISTA
				($s=1$)		($s=2$)		($s=3$)		($s=4$)		($s=5$)
time		$6.8225$		$0.5255$		$0.7784$		$1.0089$		$1.4913$		$1.4991$
iter		$50$		$50$		$50$		$50$		$50$		$50$
$\eta$		$0.2184$		$0.2649$		$0.2329$		$0.2212$		$0.2186$		$0.2182$
$\gamma$		$0.0115$		$0.0180$		$0.0123$		$0.0116$		$0.0115$		$0.0115$

Table 2. Table 2 : Numerical results for FISTA and sFISTA for ‘satellite’

		FISTA		sFISTA		sFISTA		sFISTA		sFISTA		sFISTA
				($s=1$)		($s=2$)		($s=3$)		($s=4$)		($s=5$)
time		$6.8618$		$0.5593$		$0.7870$		$0.9975$		$1.2309$		$1.6005$
iter		$50$		$50$		$50$		$50$		$50$		$50$
$\eta$		$0.2740$		$0.3551$		$0.2979$		$0.2783$		$0.2757$		$0.2744$
$\gamma$		$0.0129$		$0.0256$		$0.0152$		$0.0129$		$0.0130$		$0.0129$

Table 3. Table 3 : Numerical results for FISTA and sFISTA for ‘pattern1’

		FISTA		sFISTA		sFISTA		sFISTA		sFISTA		sFISTA
				($s=1$)		($s=2$)		($s=3$)		($s=4$)		($s=5$)
time		$7.5562$		$0.6701$		$0.8669$		$1.0097$		$1.2027$		$1.5385$
iter		$50$		$50$		$50$		$50$		$50$		$50$
$\eta$		$0.0607$		$0.4781$		$0.2171$		$0.1417$		$0.0813$		$0.0689$
$\gamma$		$0.0087$		$0.0336$		$0.0137$		$0.0135$		$0.0108$		$0.0089$

Table 4. Table 4 : Numerical results for FISTA and sFISTA for ‘ppower’

		FISTA		sFISTA		sFISTA		sFISTA		sFISTA		sFISTA
				($s=1$)		($s=2$)		($s=3$)		($s=4$)		($s=5$)
time		$7.1086$		$0.5863$		$0.8070$		$1.0227$		$1.2676$		$1.4446$
iter		$50$		$50$		$50$		$50$		$50$		$50$
$\eta$		$0.0968$		$0.2725$		$0.1478$		$0.1160$		$0.1097$		$0.0985$
$\gamma$		$0.0096$		$0.0226$		$0.0207$		$0.0139$		$0.0132$		$0.0096$

Table 5. Table 5 : Values of t ​ r ​ a ​ t ​ i ​ o 𝑡 𝑟 𝑎 𝑡 𝑖 𝑜 tratio for PRblurdefocus with different blurring levels and types of noise

tratio	case1	case2	case3	case4	case5	case6	case7	case8	case9
$\mathrm{Blurlevel}$	‘mild’	‘mild’	‘mild’	‘medium’	‘medium’	‘medium’	‘severe’	‘severe’	‘severe’
Noise type	‘gauss’	‘laplace’	‘multi’	‘gauss’	‘laplace’	‘multi’	‘gauss’	‘laplace’	‘multi’
‘pattern1’	$0.1944$	$0.2016$	$0.1995$	$0.2371$	$0.2019$	$0.2052$	$0.2064$	$0.2027$	$0.1947$
‘pattern2’	$0.1998$	$0.1971$	$0.1932$	$0.1982$	$0.2032$	$0.2007$	$0.2076$	$0.1979$	$0.1994$
‘ppower’	$0.2021$	$0.2094$	$0.2373$	$0.1950$	$0.2009$	$0.1943$	$0.1994$	$0.1957$	$0.2022$
‘smooth’	$0.2014$	$0.1983$	$0.1995$	$0.2068$	$0.1970$	$0.1999$	$0.1971$	$0.1965$	$0.2014$
‘dot2’	$0.1968$	$0.1993$	$0.2171$	$0.2011$	$0.2024$	$0.1986$	$0.2028$	$0.2028$	$0.2074$
‘dotk’	$0.2041$	$0.1995$	$0.1982$	$0.2021$	$0.2025$	$0.2143$	$0.1987$	$0.2016$	$0.2101$
‘satellite’	$0.1987$	$0.1986$	$0.1989$	$0.2164$	$0.2088$	$0.2249$	$0.2046$	$0.1996$	$0.2049$
‘hst’	$0.2006$	$0.1992$	$0.2027$	$0.1968$	$0.2016$	$0.2064$	$0.2165$	$0.2014$	$0.2034$

Table 6. Table 6 : Values of t ​ r ​ a ​ t ​ i ​ o 𝑡 𝑟 𝑎 𝑡 𝑖 𝑜 tratio for PRblurshake with different blurring levels and types of noise

tratio	case1	case2	case3	case4	case5	case6	case7	case8	case9
$\mathrm{Blurlevel}$	‘mild’	‘mild’	‘mild’	‘medium’	‘medium’	‘medium’	‘severe’	‘severe’	‘severe’
Noise type	‘gauss’	‘laplace’	‘multi’	‘gauss’	‘laplace’	‘multi’	‘gauss’	‘laplace’	‘multi’
‘pattern1’	$0.1984$	$0.2024$	$0.2172$	$0.1961$	$0.1984$	$0.1958$	$0.1973$	$0.2027$	$0.2138$
‘pattern2’	$0.1995$	$0.2054$	$0.2167$	$0.1986$	$0.1974$	$0.1992$	$0.2010$	$0.1996$	$0.2219$
‘ppower’	$0.2001$	$0.2028$	$0.2155$	$0.2081$	$0.2027$	$0.2195$	$0.2083$	$0.2113$	$0.2348$
‘smooth’	$0.2025$	$0.2015$	$0.2022$	$0.2061$	$0.2083$	$0.2223$	$0.2168$	$0.2022$	$0.2169$
‘dot2’	$0.1953$	$0.2052$	$0.2107$	$0.2141$	$0.1996$	$0.2034$	$0.2244$	$0.2350$	$0.1972$
‘dotk’	$0.2114$	$0.2083$	$0.2148$	$0.2125$	$0.2051$	$0.2189$	$0.2049$	$0.2141$	$0.2081$
‘satellite’	$0.2546$	$0.2129$	$0.2167$	$0.1965$	$0.2550$	$0.2123$	$0.2118$	$0.2210$	$0.2206$
‘hst’	$0.1991$	$0.2047$	$0.2233$	$0.2071$	$0.1984$	$0.2577$	$0.2026$	$0.2025$	$0.2160$

Equations

$$b=Ax+e,$$

$${\hat{x}}_{\mathrm{LS}}=\arg \underset{x}{\min }\frac{1}{2}\Vert Ax-b{\Vert }_2^2.$$

$${\hat{x}}_{\mathrm{TIK}}=\arg \underset{x}{\min }\frac{1}{2}\Vert Ax-b{\Vert }_2^2+\frac{{\lambda }^2}{2}\Vert Rx{\Vert }_2^2.$$

$$\underset{x}{\min }\Phi (x)=\frac{1}{2}\Vert Ax-b{\Vert }_2^2+\frac{{\lambda }^2}{2}\Vert x{\Vert }_2^2$$

$$K\otimes H=\left[\begin{array}[]{ccc}k_{11}H&\cdots&k_{1n}H\\ \vdots&&\vdots\\ k_{n1}H&\cdots&k_{nn}H\end{array}\right].$$

$$A\approx A_s=\sum _{i=1}^s{K}_i\otimes H_i,$$

$$b=A_s{x}_s+e$$

$$\underset{x_s}{\min }{\Phi }_s(x)=\frac{1}{2}\Vert A_s{x}_s-b{\Vert }_2^2+\frac{{\lambda }^2}{2}\Vert x_s{\Vert }_2^2$$

$$B=\sum _{i=1}^s{H}_{i}X{K}_i^T+E.$$

$$\underset{X}{\min }\frac{1}{2}\Vert \sum _{i=1}^s{H}_{i}X{K}_i^T-B{\Vert }_F^2+\frac{{\lambda }^2}{2}\Vert X{\Vert }_F^2,$$

$$A=\sum _{i=1}^{\min (m,n)}{K}_i\otimes H_i,$$

$$A\approx A_s=\sum _{i=1}^s{K}_i\otimes H_i.$$

$${\Vert A-\sum _{i=1}^s{K}_i\otimes H_i\Vert }_F={\Vert \bar{P}-\sum _{i=1}^s{\sigma }_i{u}_i{v_i}^T\Vert }_F,$$

$$K_i=\mathrm{toep}(k_i,l)\text{and}{H}_i=\mathrm{toep}(h_i,q).$$

$$\operatorname{toep}(c,4)=\left[\begin{array}[]{ccccccccc}c_{4}&&c_{3}&&c_{2}&&c_{1}&&0\\ c_{5}&&c_{4}&&c_{3}&&c_{2}&&c_{1}\\ 0&&c_{5}&&c_{4}&&c_{3}&&c_{2}\\ 0&&0&&c_{5}&&c_{4}&&c_{3}\\ 0&&0&&0&&c_{5}&&c_{4}\\ \end{array}\right]\qquad\text{with}\qquad c=\left[\begin{array}[]{c}c_{1}\\ c_{2}\\ c_{3}\\ c_{4}\\ c_{5}\\ \end{array}\right].$$

$$K_i=\mathrm{toep}(k_i,l)+\mathrm{hank}(k_i,l)\text{and}{H}_i=\mathrm{toep}(h_i,q)+\mathrm{hank}(h_i,q),$$

$$\operatorname{hank}(c,3)=\left[\begin{array}[]{ccccccccc}c_{4}&&c_{5}&&0&&0&&0\\ c_{5}&&0&&0&&0&&0\\ 0&&0&&0&&0&&0\\ 0&&0&&0&&0&&c_{1}\\ 0&&0&&0&&c_{1}&&c_{2}\\ \end{array}\right]\qquad\text{with}\qquad c=\left[\begin{array}[]{c}c_{1}\\ c_{2}\\ c_{3}\\ c_{4}\\ c_{5}\\ \end{array}\right].$$

$$\underset{x}{\min }\{F(x)=f(x)+g(x)\},$$

$$Q_{L(f)}(x,y):=f(y)+⟨x-y,\nabla f(y)⟩+\frac{L(f)}{2}\Vert x-y{\Vert }_2^2+g(x)$$

$$x_k=p_{L_f}(y_k)$$

$$=\arg \underset{x}{\min }\left\{⟨x,A^T(A{y}_k-b)⟩+\frac{L_f}{2}\Vert x{\Vert }_2^2+L_f⟨x,y_k⟩+\frac{{\lambda }^2}{2}\Vert x{\Vert }_2^2\right\},$$

$$\displaystyle=\arg\min_{x}\ \left\{\frac{L_{f}+{\lambda}^{2}}{2}\left\Arrowvert x-\frac{1}{L_{f}+{\lambda}^{2}}\left(L_{f}y_{k}-A^{T}(Ay_{k}-b)\right)\right\Arrowvert_{2}^{2}\right\},$$

$$=\frac{1}{L_f+{\lambda }^2}\left(L_f{y}_k-A^T(A{y}_k-b)\right).$$

$$x_k=\frac{1}{L_f+{\lambda }^2}\left(L_f{y}_k-A^T(A{y}_k-b)\right).$$

$$t_{k+1}=\frac{1+\sqrt{1+4{t_k}^2}}{2}.$$

$$y_{k+1}=x_k+\frac{t_k-1}{t_{k+1}}(x_k-x_{k-1}).$$

$$F(x_k)-F(x_F^{\ast })\le \frac{2L(f)\Vert x_0-x_F^{\ast }{\Vert }_2^2}{(k+1{)}^2},$$

$$Q_L(X,Y):=$$

$$+\frac{L}{2}\Vert X-Y{\Vert }_F^2+\frac{{\lambda }^2}{2}\Vert X{\Vert }_F^2,$$

$$Y_{k+1}=X_k+\frac{t_k-1}{t_{k+1}}(X_k-X_{k-1}),$$

$$X_{k+1}$$

$$=\arg \underset{X}{\min }\left\{{⟨X,\sum _{j=1}^r{H}_j^T\left(\sum _{i=1}^r{H}_i{Y}_{k+1}{K}_i^T-B\right)K_j⟩}_F+\frac{L}{2}\Vert X{\Vert }_F^2+L⟨X,Y_{k+1}{⟩}_F+\frac{{\lambda }^2}{2}\Vert X{\Vert }_F^2\right\},$$

$$\displaystyle=\arg\min_{X}\ \left\{\frac{L+{\lambda}^{2}}{2}\left\Arrowvert X-\frac{1}{L+{\lambda}^{2}}\left(LY_{k+1}-\sum\limits_{j=1}^{r}H_{j}^{T}\left(\sum\limits_{i=1}^{r}H_{i}Y_{k+1}K_{i}^{T}-B\right)K_{j}\right)\right\Arrowvert_{F}^{2}\right\},$$

Full text

1 Introduction

Image restoration problems have a wide range of important applications, such as digital camera and video, microscopy, meidcal imaging, etc.. Image restoration is the process of reconstructing an image of an unknown scene from an observed image, where the distortion can arise from many sources, such as motion blurs, out of focus lens, or atmospheric turbulence. Suppose there is an exact image of being all black except for a single bright pixel. If we take a picture of this image, then the distortion operation will cause the single bright pixel to be spread over its neighboring pixels. This single bright pixel is called a point source, and the function that describes the distortion and the resulting image of the point source is called the point spread function (PSF) [12]. Mathematically, the distortion can be represented by a PSF. If the PSF is the same regardless of the location of the point source, it is called spatially invariant. Throughout this paper, we assume the PSF under consideration is always spatially invariant.

A spatially invariant image restoration problem can be modeled by a linear inverse problem of the following form

[TABLE]

where $A\in{\mathbb{R}}^{N\times N}$ is a blurring matrix constructed from the PSF, $e\in{\mathbb{R}}^{N}$ is a vector representing additive noise, $b\in{\mathbb{R}}^{N}$ represents the distorted image and $x\in{\mathbb{R}}^{N}$ denotes the unknown true image to be estimated. The matrix $A$ is usually very ill-conditioned in these image restoration problems.

A classical way to solve (1) is by the least squares (LS) approach[6], whose solution takes the following form

[TABLE]

However, when $A$ is ill-conditioned, the LS solution usually has a huge norm and is thus meaningless[12]. In order to compute a decent approximation to $x$, it is necessary to employ some form of regularization. The basic idea of regularization is to replace the original ill-conditioned problem with a “nearby” well-conditioned problem whose solution is close to the orignal solution. Tikhonov regularization [19] is one of the most popular regularization techniques, where a quadratic penalty is added to the object function

[TABLE]

The second term in the above equation is a regularization term, which controls the norm (or seminorm) of the solution. The reguarization parameter $\lambda\!>\!0$ controls “smoothness” of the regularized solution. Typical choices of $R$ include an identity matrix and a matrix approximating the first or second order derivative operator [11, 14, 13].

In this paper, we choose $R$ as an identity matrix and consider the following minimization model

[TABLE]

In many applications, such as image restoration, it may also be important to include convex constraints (e.g., $x\geq 0$) on the solution.

Numerous algorithms proposed in the literature can be used to solve (2) with convex constraints. One of them is the interior point method [5, 18]. However, image restoration problems often involve dense matrix data, which will hamper the effectiveness of the interior point method. Another popular class of methods for solving (2) are gradient-based algorithms [2, 3, 23]. Although these algorithms are relatively inexpensive at each iteration, they often suffer from slow convergence. One recent development is the fast iterative shrinkage-thresholding algorithm (FISTA) [1], which was proposed to solve nonsmooth convex optimization problems. FISTA preserves the computational simplicity and has a fast global convergence rate. Thus, FISTA becomes quite attractive for solving large-scale problems. Although problem (2) does not involve any nonsmooth term, incorporating convex constraints is important in image deblurring applications. Moreover, in some situations $l_{1}$-based regularization has to be exploited to enforce sparsity in the solution. We plan to apply the proposed method to solve this class of nonsmooth optimization problems in the future. In this paper, we will first fully take advantage of the hidden structures of the blurring matrix $A$ and improve the efficiency of the FISTA framework for solving the smooth optimization problem (2).

Since the blurring model is essentially a convolution, the first structure to be exploited is the Kronecker product structure. Assume $K\in{\mathbb{R}}^{n\times n}$ and $H\in{\mathbb{R}}^{m\times m}$, the Kronecker product of these two matrices is defined as

[TABLE]

For the blurring operator $A$ in (1), it has been shown that $A$ can be approximated by a matrix $A_{s}$ as follows [15, 17]

[TABLE]

where $K_{i}\in{\mathbb{R}}^{n\times n}$, $H_{i}\in{\mathbb{R}}^{m\times m}$ with $N=mn$. The error between the blurring matrix and the Kronecker product approximation can be easily controlled. In addition, these $K_{i}$ and $H_{i}$ are not general dense matrices but structured matrices (Toeplitz, Hankel, etc. [20, 21, 22]). We will give more details on the error between $A$ and $A_{s}$ and the structures of $K_{i}$ and $H_{i}$ in Section 2.

Consequently, the solution of (1) can be approximated by the following problem

[TABLE]

and equivalently the solution of (2) can be approximated by solving the optimization problem

[TABLE]

From the numerical examples in Section 4, we can see that $x_{s}$ from (6) and $x$ from (2) can provide indistinguishable image restoration results. This is because the original ill-posed problem (2) only requires a numerical solution $x$ with relatively low accuracy. As long as the difference between $A_{s}$ and $A$ falls below a certain level, which can be easily met with only a small value of $s$ in (4), $x_{s}$ from (6) and $x$ from (2) can reach the same level of accuracy. This phenomenon is analyzed in Theorem 4 in Section 3 and verified by the numerical experiments in Section 4.

If $b=\mathrm{vec}(B),\ x_{s}=\mathrm{vec}(X)\ \text{and}\ e=\mathrm{vec}(E)$, where $\mathrm{vec}(X)$ represents a column vector obtained from vectorizing a matrix $X$ (i.e. columns of $X$ are stacked one after the other), then (5) can be rewritten equivalently as

[TABLE]

It is straightforward to derive the corresponding Tikhonov regularized minimization model as follows

[TABLE]

where $\|\cdot\|_{F}$ denotes the Frobenius norm. (8) has several advantages over the original optimization problem (2). First, (8) benefits from the Kronecker product structure of $A_{s}$ and can exploit more computationally intensive matrix-matrix operations. In addition, all the matrices $H_{i}$ and $K_{i}$ are structured matrices, which enables fast matrix-vector multiplications at each iteration. Second, the summation of $s$ terms in (8) can be performed independently and enables (8) to reach superior parallel efficiency when implemented on modern high performance computing architectures. Some work has been done to exploit matrix equation structures for iterative methods to solve inverse problems of the form (7); see, for example, [4, 7, 8, 24]. In this paper, we propose the structured FISTA (sFISTA) method. It gains its efficiency by exploiting both the Kronecker product structure of $A$ as well as the structures from $K_{i}$ and $H_{i}$. The convergence rate of sFISTA can be of the same order as FISTA under mild conditions.

The remaining sections are organized as follows. In Section 2, we describe how to approximate the blurring matrix $A$ into the sum of a few of Kronecker products. In Section 3, we first briefly review the FISTA framework and then propose the sFISTA method. We also show that sFISTA for (8) is equivalent to FISTA for (6) and derive the convergence and complexity analysis of sFISTA for (8). Some numerical examples are provided in Section 4 and the concluding remarks are drawn in Section 5.

2 Kronecker Decomposition

Consider a 2-D spatially-invariant image restoration problem. It was shown in [16] that three different structures of the blurring matrix $A$ commonly occur. If the zero boundary condition (corresponding to assuming the values of $x$ outside the domain of consideration are zero) is applied, $A$ will be a block-Toeplitz-Toeplitz-block (BTTB) matrix. On the other hand, if the periodic boundary condition (corresponding to the case that the image outside the domain of consideration is a repeat of the image inside in all directions) is used, $A$ becomes a block-circulant-circulant-block (BCCB) matrix. Finally, $A$ would be block-Toeplitz-plus-Hankel with Toeplitz-plus-Hankel-blocks (BTHTHB) if the reflective boundary condition (corresponding to a reflection of the original scene at the boundary) is utilized. In any case, the matrix $A$ can always be approximated as the sum of a few Kronecker products. Since the periodic boundary condition often cause severe ringing artifacts near image borders, only the other two cases are considered in the remaining sections.

In practice, the PSF for images with $m\times n$ pixels is often stored as an $m\times n$ array $P$. When $P$ represents the image of a single bright pixel, the process of taking a picture of such an image is equivalent to computing one column of matrix $A$ with column index $t$, where $t$ depends on the location of the point source. Thus, the structure of $A$ is completely determined by that of $P$. More specifically, suppose $P$ has the SVD decomposition $P=U\Sigma V^{T}$. Let $u_{i}$ and $v_{i}$ be the $i$th columns of matrices $U$ and $V$, respectively and $\sigma_{1}\geq\sigma_{2}\geq\cdots\geq\sigma_{\min{(m,n)}}$ be the singular values of $P$. It has been shown that $A$ then admits the following Kronecker decomposition [15, 17]

[TABLE]

where $K_{i}$ and $H_{i}$ are matrices defined based on $u_{i}$, $v_{i}$, $\sigma_{i}$ and boundary conditions. More details on the structure of $K_{i}$ and $H_{i}$ will be provided at the end of this section. Because the singular values of $P$ decay quickly in realistic applications, (9) can be further truncated by keeping only the first $s$ terms

[TABLE]

The approximation error introduced in (10) has been well studied in [15, 17]. The analysis in [15, 17] shows that the distance between $A$ and $A_{s}$ is related to the approximation error of a truncated SVD decomposition of a matrix $\bar{P}$, which is summarized in the following theorem for the square PSF case.

Theorem 1**.**

[17, Theorem 3.1]**. Assume the blurring matrix $A$ is constructed from a PSF $P$ with center $p_{lq}$ located at $(l,q)$, then for both zero boundary condition and reflective boundary condition, we have

[TABLE]

*where $\bar{P}=W_{a}PW_{b}$ with $W_{a}=\mathrm{diag}\left([\ \sqrt{n-l+1}\ \cdots\ \sqrt{n-1}\ \sqrt{n}\ \sqrt{n-1}\ \cdots\ \sqrt{l}\ ]^{T}\right)$, $W_{b}=\mathrm{diag}$

$\left([\ \sqrt{n-q+1},\ \cdots\ \sqrt{n-1}\ \sqrt{n}\ \sqrt{n-1}\ \cdots\ \sqrt{q}\ ]^{T}\right)$ for the zero boundary condition case and $\bar{P}=RPR^{T}$ with $R$ is the Cholesky factor of the symmetric Toeplitz matrix with its first row as $[n,1,0,1,0,1,\cdots]$ for the reflective boundary condition case. Here $\sum\limits_{i=1}^{s}\sigma_{i}u_{i}{v_{i}}^{T}$ is the summation of the first $s$ terms in the SVD decomposition of $\bar{P}$.*

Since the singular values of $\bar{P}$ (as well as $P$) decay quickly to zero for most PSFs, Theorem 1 guarantees that even a small $s$ in (10) could lead to very accurate approximation. Numerical experiments in Section 4 show that taking $s$ as small as $5$ is enough for the image restoration applications under consideration.

At the end of this section, let us take a look at the structure of $K_{i}$ and $H_{i}$. If the zero boundary condition is used, $K_{i}$ and $H_{i}$ have the following Toeplitz structure

[TABLE]

In the above equations, $k_{i}=(\sqrt{\sigma_{i}}u_{i})./\operatorname{diag}{(W_{a})}$, $h_{i}=(\sqrt{\sigma_{i}}v_{i})./\operatorname{diag}{(W_{b})}$, where $./$ denotes point-wise division. And $\operatorname{toep}(c,j)$ denotes a banded Toeplitz matrix whose $j$th column is equal to $c$. For example

[TABLE]

On the other hand, if the reflective boundary condition is applied, $K_{i}$ and $H_{i}$ are equal to the linear combinations of a Toeplitz matrix and a Hankel matrix:

[TABLE]

where $k_{i}=\sqrt{\sigma_{i}}R^{-1}u_{i}$, $h_{i}=\sqrt{\sigma_{i}}R^{-1}v_{i}$ and $\operatorname{hank}(c,j)$ denotes a banded Hankel matrix whose first row and last column are defined by $[c_{j+1},\cdots,c_{n},0,\cdots,0]$ and $[0,\cdots,0,c_{1},\cdots,c_{j-1}]^{T}$, respectively. For example

[TABLE]

The sFISTA to be introduced in the next section will benefit from the fast Toeplitz/Hankel matrix-vector product algorithms when multiplying $K_{i}$ and $H_{i}$ with vectors at each iteration.

3 Structured FISTA

In this section, we will first review the FISTA framework for solving (2) and then propose sFISTA for solving (8). We can prove that the proposed sFISTA for solving (8) is equivalent to FISTA for solving (6). A detailed error analysis has also been conducted to show that the computational accuracy of sFISTA can reach the same level as that of FISTA under mild conditions. Finally, we compare the computational complexity of sFISTA for solving (8) and FISTA for solving (2) and show that sFISTA is more efficient in both serial and parallel computing environments.

3.1 FISTA: A fast iterative shrinkage-thresholding algorithm

A fast iterative shrinkage-thresholding algorithm (FISTA) was first proposed in [1] to solve the following general nonsmooth convex optimization model

[TABLE]

where $f:{\mathbb{R}^{N}\rightarrow\mathbb{R}}$ is a smooth convex function of the type $\mathrm{C^{1,1}}$ and $g:{\mathbb{R}^{N}\rightarrow\mathbb{R}}$ is a continuous convex function which is possibly nonsmooth. The basic idea of FISTA is that at each iteration, after getting the current iteration point $x_{k}$, an additional point $y_{k+1}$ is chosen as the linear combination of the current iteration point $x_{k}$ and the previous iteration point $x_{k-1}$. The next iteration point $x_{k+1}$ is then set as the unique minimizer $p_{L(f)}(y_{k+1})$ of the quadratic approximation $Q_{L(f)}(x,y_{k+1})$ of $F(x)$ at $y_{k+1}$ with

[TABLE]

and $L(f)$ being the Lipschitz constant of $\nabla f$. For more details about FISTA, one can refer to [1].

Obviously, (2) is a special instance of problem (16) if we let $f(x)=\frac{1}{2}\|Ax-b\|_{2}^{2}$ and $g(x)=\frac{{\lambda}^{2}}{2}\|x\|_{2}^{2}$. In this case, the (smallest) Lipschitz constant of the gradient $\nabla f$ is $L_{f}=\lambda_{\mathrm{max}}(A^{T}A)$. Simple calculations lead to

[TABLE]

See Algorithm 1 for a description of FISTA for (2).

As can be seen from Algorithm 1, the total computational cost of FISTA is dominated by matrix-vector multiplications associated with $A$ and $A^{T}$ at Step $1$. Other steps only involve inexpensive vector and scalar operators. Despite its simplicity, FISTA enjoys a fast global convergence rate, which is summarized in Theorem 2.

Theorem 2**.**

[1, Theorem 4.4]* Let $\{x_{k}\}$, $\{y_{k}\}$ be generated by FISTA. Then for any $k\geq 1$*

[TABLE]

where $x_{F}^{*}$ is the solution of (16).

It is well known that many first order algorithms have very slow convergence rate. From Theorem 2, we can see that FISTA is different from classical first order methods in the sense that it preserves a fast global convergence rate $O(1/k^{2})$. That is, in order to obtain a numerical solution $x$ such that $F(x)-F(x_{F}^{*})\leq\epsilon$, the number of iterations required by FISTA is at most $\frac{\sqrt{2L(f)}\|x_{0}-x_{F}^{*}\|_{2}}{\sqrt{\epsilon}}-1$. In the next section, we will propose the sFISTA which is more efficient for solving (8).

3.2 Accelerating FISTA by exploiting structures

In this section, we will show how to adapt the FISTA framework to solve (8) by exploiting the two hidden structures. We first use a Kronecker product approximation $A_{s}$ of the coefficient matrix $A$ to introduce problem (6), which can be equivalently transformed into a matrix problem (8). Consider the following quadratic approximation of the objective function of (8) at a given point $Y$:

[TABLE]

where $L=\lambda_{\mathrm{max}}(A^{T}A)$ is the Lipschitz constant of the gradient of the first term in the object function of (8). Similar to FISTA, we choose the unique minimizer of the quadratic approximation at point $Y_{k+1}$, which is the linear combination of $X_{k}$ and $X_{k-1}$, as the new iteration point $X_{k+1}$. Mathematically, we set

[TABLE]

where $t_{k}$ and $t_{k+1}$ are parameters updated in the same way as FISTA to make sFISTA maintain the same convergence rate as FISTA for solving (8) and compute $X_{k+1}$ as

[TABLE]

Basic steps of sFISTA for (8) are summarized in Algorithm 2.

Compared with Algorithm 1, there are several major differences between sFISTA and FISTA. First of all, the computational cost of Algorithm 1 is dominated by matrix-vector multiplications while Algorithm 2 can benefit from more computationally intensive matrix-matrix multiplications. Moreover, since $H_{i}$ and $K_{i}$ are all structured matrices (Toeplitz, Hankel, etc.), we can further exploit their fast matrix-vector multiplications at Step $1$ in Algorithm 2. Second, Algorithm 2 decomposes the computation of $X_{k}$ as the summation of $s$ terms, which can be computed independently. Therefore, we can easily explore two levels of parallelism at each iteration in Algorithm 2. The first level corresponds to the structured matrix-vector multiplications with multiple vectors and the second level comes from the summation of $s$ terms. This property enables Algorithm 2 to reach superior parallel performance when implemented on modern high performance architectures. Finally, we can prove that sFISTA for (8) is equivalent to FISTA for (6), which guarantees the fast convergence.

Theorem 3**.**

sFISTA for (8) and FISTA for (6) provide the same output as long as their initial points satisfy $x_{0}=\mathrm{vec}(X_{0})$. Mathematically, suppose $\{X_{k}\}$, $\{Y_{k}\}$ are generated by sFISTA for (8) and $\{x_{k}\}$, $\{y_{k}\}$ are obtained by FISTA for (6), then we have $x_{k}=\mathrm{vec}(X_{k})$ and $y_{k}=\mathrm{vec}(Y_{k})$.

Proof.

To prove the desired results, we first review two important properties of Kronecker products, which will be used in the proof below.

[TABLE]

where $H$, $K$ and $Z$ are matrices of appropriate dimensions. Recall the frameworks of two algorithms, to prove they provide the same output, we only have to show that both algorithms are equivalent at Step $1$. Specifically, we just have to prove

[TABLE]

Utilizing the Kronecker product properties mentioned above, we can get

[TABLE]

from which we can derive that $x_{k}=\mathrm{vec}(X_{k})$ and $y_{k}=\mathrm{vec}(Y_{k})$. ∎

It is worth pointing out that sFISTA for (8) is only equivalent to FISTA for (6) due to the Kronecker product approximation error. The total computational error of sFISTA for solving (8) comes from two places: the Kronecker product approximation to $A$ and the iterative procedure of sFISTA. The following theorem analyzes the effect of these two kinds of errors on the accuracy of the final computed result.

Theorem 4**.**

Assume $x^{*}$ and $x_{s}^{*}$ are the exact solutions of (2) and (6) respectively, $\{{X}_{k},{Y}_{k}\}$ is the sequence obtained by sFISTA for (8). Denote $\tilde{x}_{k}=\mathrm{vec}(X_{k})$. If the singular values of $A^{T}A_{s}+\lambda^{2}I$ have a lower bound and the Kronecker product approximation $A_{s}$ satisfies $\|A-A_{s}\|_{F}=\epsilon_{s}$, then we have for any $k\geq 1$

[TABLE]

where $c_{0}$ is a positive constant independent of $k$.

Proof.

Since $x^{*}$ and $x_{s}^{*}$ are the exact solutions of (2) and (6), respectively, from their optimality conditions we have

[TABLE]

which implies that

[TABLE]

It is easy to see

[TABLE]

For the first term we have

[TABLE]

where we use the fact that $\|A\|_{2}$, $\|A_{s}\|_{2}$, $\|\tilde{x}_{k}\|_{2}$ are bounded.

From Theorem 3 we know that sFISTA for (8) is equivalent to FISTA for (6), which implies that the second term satisfies

[TABLE]

To estimate the last term, we first prove the following fact. From (22) we get

[TABLE]

which implies that

[TABLE]

Then we have

[TABLE]

where we utilize the boundedness of $\|A\|_{2}$, $\|A_{s}\|_{2}$, $\|b\|_{2}$, $\|x^{*}\|_{2}$, $\|x_{s}^{*}\|_{2}$ and the assumption that the singular values of $A^{T}A_{s}+\lambda^{2}I$ have a lower bound.

Then for the last term we have

[TABLE]

where we utilize the boundedness of $\|A\|_{2}$, $\|b\|_{2}$, $\|x_{s}^{*}\|_{2}$ and the fact that $\|x_{s}^{*}-x^{*}\|_{2}\leq\tilde{c}\epsilon_{s}$.

Based on the above analysis above for the three terms in (24), it follows that

[TABLE]

Let $c_{0}=c_{1}+c_{2}$, then the desired result (21) follows. ∎

Theorem 4 shows that the error $|\Phi(\tilde{x}_{k})-\Phi(x^{*})|$ from sFISTA is bounded by two terms: $\frac{2L\|\tilde{x}_{0}-x_{s}^{*}\|_{2}^{2}}{(k+1)^{2}}$ and $c_{0}\epsilon_{s}$. The first term decreases as the iteration proceeds while the second term remains constant during the iteration. In order to let the total error $|\Phi(\tilde{x}_{k})-\Phi(x^{*})|$ fall below a threshold $\epsilon$, we need to make both terms smaller than $\epsilon$. As discussed before, since only a relatively large $\epsilon$ is necessary in these ill-posed inverse problems, a small $s$ would be enough to guarantee $c_{0}\epsilon_{s}<\epsilon$. In this sense, the convergence of sFISTA is dominated by the first term $\frac{2L\|\tilde{x}_{0}-x_{s}^{*}\|_{2}^{2}}{(k+1)^{2}}$ and behaves in a similar way as FISTA.

As an example, we plot the singular values of the matrices $P$ and $\bar{P}$ from the test image ‘hst’ (See Example 1 in Section 4 for more details about this image) in Figure 1. It is easy to see that the singular values of both matrices decay quickly to zero. For example, the ratio of the sixth largest singular value of $P$ to the largest one is only $6.47e\!-\!2$ and the ratio of the tenth largest singular value of $P$ to the largest one reduces to $4.27e\!-\!2$. These patterns can also be observed in other test examples.

3.3 Complexity Analysis

In this section, we consider the computational complexity of sFISTA for (8) (Algorithm 2) and FISTA for (2) (Algorithm 1). If we ignore the structures in $A\in\mathcal{R}^{mn\times mn}$, $K_{i}\in\mathcal{R}^{n\times n}$, $H_{i}\in\mathcal{R}^{m\times m}$ and assume that they are all general dense matrices, then the cost of Step $1$ in Algorithm 1 and Algorithm 2 would be $O\left(m^{2}n^{2}\right)$ and $O\left(s(nm^{2}+mn^{2})\right)$, respectively. When $s$ is much smaller than $m$ and $n$, which is the case for the applications under consideration in this paper, Algorithm 2 is definitely faster than Algorithm 1.

Recall that the blurring matrix $A$ and matrices $K_{i}$ and $H_{i}$ from the Kronecker product approximation of $A$ all have specific structures. As the matrix size becomes big enough, these structures will enable us to use fast Fourier transforms (FFTs) to accelerate matrix-vector multiplications encountered in both algorithms. For example, when zero boundary condition is used, $A$ is a block-Toeplitz-Toeplitz-block (BTTB) matrix and $K_{i}$, $H_{i}$ are Toeplitz matrices. In this case, the matrix-vector multiplication at step $1$ in Algorithm 1 can be performed in $O(mn\log(mn))$ with 2D FFTs, while Step $1$ in Algorithm 2 can be done with 1D FFTs in $O\left(smn\log(mn)\right)$. When reflective boundary condition is utilized, $A$ is a block-Toeplitz-plus-Hankel with Toeplitz-plus-Hankel-blocks (BTHTHB) matrix and $K_{i}$, $H_{i}$ can be represented as the sum of a Toeplitz matrix and a Hankel matrix. In this case, the computational complexities of Step $1$ in both algorithms are still of the same order as in the zero boundary condition case. Although Algorithm 2 has the same complexity as Algorithm 1, it is important to notice that Algorithm 2 is actually much more attractive when implemented on high performance architectures for a number of reasons. First of all, as discussed in the previous section, Algorithm 2 can easily exploit two levels of parallelism, which is crucial for fully taking advantage of the multilelvel parallelism offered by the current architectures. Second, parallel 1D FFTs are known to scale better than parallel 2D FFTs. Thus, Algorithm 2 is more computationally efficient than Algorithm 1 for solving large scale problems.

4 Numerical Results

In this section, we provide some numerical examples to demonstrate the performance of sFISTA for solving (8). All the algorithms were implemented with MATLAB and the experiments were performed on a Macbook Air with Intel Core i7 CPU (2.2 GHz). The following notations will be used throughout the section:

• $\bullet$

$s$: the number of terms in the Kronecker product approximation;

• $\bullet$

$b$: the data vector;

• $\bullet$

$\mathrm{noise}$: the vector of perturbations;

• $\bullet$

$b_{n}$: the noisy data $b_{n}=b+\mathrm{noise}$;

• $\bullet$

$\mathrm{NoiseLevel}$: relative level of noise defined as $\|\mathrm{noise}\|_{2}/\|b\|_{2}$

• $\bullet$

$\mathrm{BlurLevel}$: an indicator used to set the severity of the blur to one of the following: ‘mild’, ‘medium’ and ‘severe’;

• $\bullet$

$\eta$: the relative error $\|x-x^{*}\|/\|x^{*}\|$;

• $\bullet$

$\gamma$: the relative residual $\|r\|/\|b\|$, where $r=Ax-b$;

• $\bullet$

$iter$: the iteration number of one algorithm;

• $\bullet$

$t(\mathrm{FISTA})$ and $t(\mathrm{sFISTA})$: the CPU time (seconds) of FISTA and sFISTA, respectively;

• $\bullet$

$tratio$: an indicator defined as $tratio=\frac{t(\mathrm{sFISTA})}{t(\mathrm{FISTA})}$ to compare the efficiency of FISTA and sFISTA.

Example 1**.**

In this example, four $256\times 256$ simple test images were extracted based on functions PRblurdefocus and PRblurshake from the regularization toolbox [10]. The four test images in this example are represented by ‘hst’ (image of the Hubble space telescope), ‘satellite’ (satellite test image), ‘pattern1’ (geometrical image) and ‘ppower’ (random image with patterns of nonzero pixels) respectively, which used reflective (Neumann) boundary conditions [12]. PRblurdefocus and PRblurshake are functions simulating a spatially invariant, out-of-focus blur and spatially invariant motion blur caused by shaking of a camera, respectively. The $\mathrm{BlurLevel}$ was set to be ‘medium’ in these four tests. In addition, function PRnoise was used to add Gaussian noise with $\mathrm{NoiseLevel}=0.01$ in this example. The regularization parameters were chosen automatically by IRhybrid$\_$lsqr from [10], which is based on the hybrid bidiagonalization method presented in [9]. The Lipschitz constant was computed as an estimation of the 2-norm of the matrix $A$, which was realized by a few iterations of Lanczos bidiagonalization as implemented in HyBR [10].

We then tested FISTA for (2) and sFISTA for (8) on these four images. To show how the number of terms in the Kronecker product approximation affects the performance of sFISTA, $s$ was set to range from $1$ to $5$ in these four tests. The maximum iteration number for both algorithms was fixed at $50$. To compare the performance of FISTA and sFISTA, we report the CPU time (seconds), the relative error $\eta$ and the relative residual $\gamma$ returned by both algorithms. Their values on these four tests are tabulated in Tables 1–4.

As can be seen from Tables 1–4, sFISTA is much faster than FISTA in all test problems. As $s$ increases from $1$ to $5$, the computational time of sFISTA increases monotonically while both the relative error $\eta$ and the relative residual $\gamma$ keep decreasing. When $s$ reaches $5$, the errors of sFISTA are close enough to those of FISTA and sFISTA is still about $5$ times faster than FISTA. We would like to emphasize that sFISTA was only implemented as a serial code and we expect to see a larger speedup with a parallel implementation in the future.

We also plot the four images obtained by sFISTA in Figures 2-5. As a comparison, the true, blurred and noisy image and the image obtained by FISTA are also provided. It is easy to see that the images obtained by FISTA and sFISTA ($s=5$) seem very similar to each other.

Example 2**.**

In this set of test problems, we compare the performance of sFISTA and FISTA on eight images with different blur levels and noise levels. These problems are all extracted from functions PRblurdefocus and PRblurshake. The eight test images are represented by ‘pattern1’ (geometrical image), ‘pattern2’ (geometrical image), ‘ppower’ (random image with patterns of nonzero pixels), ‘smooth’ (very smooth image), ‘dot2’ (two small Gaussian shaped dots), ‘dotk’ ($N/2$ small Gaussian shaped dots), ‘satellite’ (satellite test image) and ‘hst’ (image of the Hubble space telescope), respectively. Each test image in this example undergoes blurring and noise-adding procedure with three different blurring levels: ‘mild’, ‘medium’, ‘severe’ and three different kinds of noise ‘gauss’ (Gaussian white noise), ‘laplace’ (Laplacian noise), ‘multiplicative’ (specific type of multiplicative noise, which is often encountered in radar and ultrasound imaging). Since we have seen from Example 1 that $s=5$ is a good choice for sFISTA, $s$ is fixed to be $5$ in this example.

We then tested sFISTA and FISTA on these eight images to compare their performance. The computational results are tabulated in Table 5 and Table 6, respectively. In both tables, there are $9$ cases corresponding to three different $\mathrm{Blurlevel}$ and three different types of noise. Since the indicator $tratio$ measures the ratio of the computational time of sFISTA to that of FISTA for solving the same problem to the same accuracy. Therefore, the smaller $tratio$ is, the more efficient sFISTA is than FISTA.

The results in Tables 5 and 6 show that sFISTA is more efficient than FISTA in all test problems. For example, when we focus on one row or one column of Table 5 or Table 6, it is easy to see that all the values of $tratio$ are less than $1$. Moreover, $tratio$ in both tables are quite close to $0.2$, which indicates that the efficiency of sFISTA does not depend on the blurring type, the test images, the blurring levels or the noise types. These results further indicate that sFISTA is not only fast but also very robust for solving the image restoration problems considered in this paper.

5 Conclusion

In this paper, we propose the structured fast iterative shrinkage-thresholding algorithm (sFISTA) for solving large scale ill-posed linear inverse problems arising from image restoration. By exploiting both the Kronecker product structure of the coefficient matrix and the pattern structure of the matrices in the Kronecker product approximation, sFISTA can significantly accelerate the computation compared to FISTA. A theoretical error analysis has been conducted to show that sFISTA can reach the same level of computational accuracy as FISTA under certain conditions. Finally, the efficiency of sFISTA is demonstrated with both a theoretical computational complexity analysis and various numerical examples coming from different applications.

The proposed sFISTA framework provide the possibility of developing new solvers for other imaging deburring problems. For example, it is possible to adapt sFISTA to solve nonsmooth optimization problems with sparsity constraints such as those using $l_{1}$-based regularization. We also plan to exploit preconditioning techniques to further reduce the iteration number and iteration time in our future work.