A tensor-based dictionary learning approach to tomographic image reconstruction

Abstract

We consider tomographic reconstruction using priors in the form of a dictionary learned from training images. The reconstruction has two stages: first we construct a tensor dictionary prior from our training data, and then we pose the reconstruction problem in terms of recovering the expansion coefficients in that dictionary. Our approach differs from past approaches in that (a) we use a third-order tensor representation for our images and (b) we recast the reconstruction problem using the tensor formulation. The dictionary learning problem is presented as a non-negative tensor factorization problem with sparsity constraints. The reconstruction problem is formulated in a convex optimization framework by looking for a solution with a sparse representation in the tensor dictionary. Numerical results show that our tensor formulation leads to very sparse representations of both the training images and the reconstructions due to the ability of representing repeated features compactly in the dictionary.

Appendix: Reconstruction solution via TFOCS

The reconstruction problem (15) is a convex, but $\Vert {\mathcal {C}} \Vert _{\mathrm {sum}}$ and $\Vert C\Vert _*$ are not differentiable which rules out conventional smooth optimization techniques. The TFOCS software [4] provides a general framework for solving convex optimization problems, and the core of the method computes the solution to a standard problem of the form

$$\begin{aligned} \text{ minimize } \quad l(A(x)-b)+h(x) , \end{aligned}$$

(17)

where the functions l and h are convex, A is a linear operator, and b is a vector; moreover l is smooth and h is non-smooth.

To solve problem (15) by TFOCS, it is reformulated as a constrained linear least squares problem:

$$\begin{aligned} \min _{{\mathcal {C}}} \frac{1}{2} \left\| \begin{pmatrix} 1/\sqrt{m} \, A \\ \delta /c \, L \end{pmatrix} \Pi \mathtt {vec}({\mathcal {D}}*{\mathcal {C}}) - \begin{pmatrix} b \\ 0 \end{pmatrix} \right\| _2^2 + \mu \, \varphi _{\nu }({\mathcal {C}}) \qquad \text {s.t.} \qquad {\mathcal {C}} \ge 0, \end{aligned}$$

(18)

where $c = \sqrt{2( M(M/p-1)+N(N/r-1) )}$. Referring to (17), $l(\cdot )$ is the squared 2-norm residual and $h(\cdot ) = \mu \, \varphi _{\nu }(\cdot )$.

The methods used in TFOCS require computation of the proximity operators of the non-smooth function h. The proximity operator of a convex function is a natural extension of the notion of a projection operator onto a convex set [12].

Let $f = \Vert {\mathcal {C}}\Vert _{\mathrm {sum}} = \Vert C\Vert _{\mathrm {sum}}$ and $g=\Vert C\Vert _*$ be defined on the set of real-valued matrices and note that ${\mathrm {dom}} f \cap {\mathrm {dom}}\, g \ne \emptyset$. For $Z \in {\mathbb {R}}^{m\times n}$ consider the minimization problem

$$\begin{aligned} \text{ minimize }_X \ f(X)+g(X) + 1/2 \Vert X-Z\Vert _{\mathrm {F}}^2 \end{aligned}$$

(19)

whose unique solution is $X = {\mathrm {prox}}_{f+g}(Z)$. While the prox operators for $\Vert C\Vert _{\mathrm {sum}}$ and $\Vert C\Vert _*$ are easily computed, the prox operator of the sum of two functions is intractable. Although the TFOCS library includes implementations of a variety of prox operators—including norms and indicator functions of many common convex sets—implementation of prox operators of the form ${\mathrm {prox}}_{f+g}(\cdot )$ is left out. Hence we compute the prox operator for $\Vert \cdot \Vert _{\mathrm {sum}}+\Vert \cdot \Vert _*$ iteratively using a Dykstra-like proximal algorithm [12], where prox operators of $\Vert \cdot \Vert _{\mathrm {sum}}$ and $\Vert \cdot \Vert _*$ are consecutively computed in an iterative scheme.

Let $\tau =\mu /q \ge 0$. For $f(X)=\tau \Vert X\Vert _{\mathrm {sum}}$ and $X \ge 0$, ${\mathrm {prox}}_f$ is the one-sided elementwise shrinkage operator

$$\begin{aligned} {\mathrm {prox}}_f(X)_{i,j} = {\left\{ \begin{array}{ll} 0 , &{}\quad X_{i,j} \ge \tau \\ X_{i,j}-\tau , &{}\quad |X_{i,j}| \le \tau \\ 0, &{}\quad X_{i,j} \le -\tau \\ \end{array}\right. } \end{aligned}$$

The proximity operator of $g(X)=\tau \Vert X\Vert _*$ has an analytical expression via the singular value shrinkage (soft threshold) operator

$$\begin{aligned} {\mathrm {prox}}_g (X) = U {\mathrm {diag}}(\sigma _i-\tau ) \, V^T, \end{aligned}$$

where $X=U \Sigma \, V^T$ is the singular value decomposition of X [9]. The computation of $\tau \Vert C \Vert _*$ can be done very efficiently since C is $sq \times r$ with $r \ll sq$.

The iterative algorithm which computes an approximate solution to ${\mathrm {prox}}_{f+g}$ is given in Algorithm 2. Every sequence $X_k$ generated by Algorithm 2 converges to the unique solution ${\mathrm {prox}}_{f+g}$ of problem (19) [12].