An algorithm for constrained one-step inversion of spectral CT data

We carry out the deriviations on $D_{\text{LSQ}}$ and $D_{\text{TPL}}$ in parallel. The local quadratic expansion for each of these data discrepancy terms about the material maps $f=f_{0}$ is

		$\displaystyle L(f)\approx L\left(f_{0}\right)+\left(f-f_{0}\right)^{\top}\nabla_{f}L\left(f_{0}\right)+\frac{1}{2}\left(f-f_{0}\right)^{\top}\nabla^{2}_{f}L\left(f_{0}\right)\left(f-f_{0}\right),$		(12)
	where
		$\displaystyle L(f)=D(c,\hat{c}(f)).$

To obtain the desired expansions, we need expressions for the gradient and Hessian of each data discrepancy. The gradient of $L_{\text{TPL}}(f)$ is derived explicitly in Appendix A; we do not show the details for the other derivations. The data discrepancy gradients are:

	$\displaystyle\nabla_{f}L_{\text{TPL}}(f)$	$\displaystyle=Z^{\top}A(f)^{\top}r(f),$		(13)
	$\displaystyle\nabla_{f}L_{\text{LSQ}}(f)$	$\displaystyle=Z^{\top}A(f)^{\top}r^{\text{(log)}}(f),$		(14)

where $r$ and $r^{\text{(log)}}$ denote the residuals in terms of counts or log counts:

	$\displaystyle r_{w\ell}(f)$	$\displaystyle=c_{w\ell}-\hat{c}_{w\ell}(f),$		(15)
	$\displaystyle r^{\text{(log)}}_{w\ell}(f)$	$\displaystyle=\log(c_{w\ell})-\log\left(\hat{c}_{w\ell}(f)\right);$		(16)

$Z$ represents the combined linear transform that accepts material maps, performs projection, and then combines the resulting sinograms to form monochromatic sinograms at energy $E_{i}$:

$$Z_{\ell i,mk}=\mu_{mi}X_{\ell k};$$

(17)

and $A(f)$ is a term that results from the gradient of the logarithm of the estimated counts $\log\hat{c}(f)$:

	$\displaystyle A_{w\ell,\ell^{\prime}i}(f)=$	$\displaystyle\frac{s_{w\ell i}\exp\left[-(Zf)_{\ell i}\right]}{\sum_{i^{\prime}}s_{w\ell i^{\prime}}\exp\left[-(Zf)_{\ell i^{\prime}}\right]}\mathbf{I}_{\ell\ell^{\prime}}$		(18)
	$\displaystyle\mathbf{I}_{\ell\ell^{\prime}}=$	$\displaystyle\begin{cases}0&\ell\neq\ell^{\prime}\\ 1&\ell=\ell^{\prime}\end{cases}.$

Using the same variable and linear transform definitions, the expressions for the two Hessians are

	$\displaystyle\nabla^{2}_{f}L_{\text{TPL}}(f)$	$\displaystyle=-Z^{\top}\operatorname{diag}\left(A(f)^{\top}r(f)\right)Z+Z^{\top}A(f)^{\top}\operatorname{diag}\left(\hat{c}(f)+r(f)\right)A(f)Z,$		(19)
	$\displaystyle\nabla^{2}_{f}L_{\text{LSQ}}(f)$	$\displaystyle=-Z^{\top}\operatorname{diag}\left(A(f)^{\top}r^{\text{(log)}}(f)\right)Z+Z^{\top}A(f)^{\top}\operatorname{diag}\left(1+r^{\text{(log)}}(f)\right)A(f)Z.$		(20)

Substituting either Eq. (19) or (20) for the Hessian and either Eq. (13) or (14) for the gradient into the Taylor expansion in Eq. (12), yields the quadratic approximation to the data discrepancy terms of interest. This quadratic is in general non-convex because both Hessian expressions can have negative values.

The key to deriving a local convex upper bound to the quadratic expansion of $L(f)$ is to split the Hessian expressions into positive and negative components. Setting the negative components to zero and substituting this thresholded Hessian into the Taylor’s expansion, yields a quadratic term with non-negative curvature. (As an aside, a tighter convex local quadratic upper bound would be attained by diagonalizing the Hessian and forming a positive semi-definite Hessian by keeping eigenvectors corresponding to only non-negative eigenvalues in the eigenvalue decomposition, but for realistic sized tomography configurations such an eigenvalue decomposition is impractical.) The algebraic steps for splitting the Hessian into positive and negative components in the form

$$\nabla^{2}_{f}L(f)=\nabla^{2}_{+}L(f)-\nabla^{2}_{-}L(f),$$

where $\nabla^{2}_{+}L(f)$ and $\nabla^{2}_{-}L(f)$ are both positive semidefinite (see Appendix B for more details). The resulting split expressions are:

	$\displaystyle\nabla^{2}_{+}L_{\text{TPL}}(f)$	$\displaystyle=Z^{\top}\operatorname{diag}\left(A(f)^{\top}r_{-}(f)\right)Z+Z^{\top}A(f)^{\top}\operatorname{diag}\left(\hat{c}(f)-r_{-}(f)\right)A(f)Z,$		(21)
	$\displaystyle\nabla^{2}_{-}L_{\text{TPL}}(f)$	$\displaystyle=Z^{\top}\operatorname{diag}\left(A(f)^{\top}r_{+}(f)\right)Z-Z^{\top}A(f)^{\top}\operatorname{diag}\left(r_{+}(f)\right)A(f)Z,$		(22)

and

	$\displaystyle\nabla^{2}_{+}L_{\text{LSQ}}(f)$	$\displaystyle=Z^{\top}\operatorname{diag}\left(A(f)^{\top}r_{-}^{\text{(log)}}(f)\right)Z+Z^{\top}A(f)^{\top}\operatorname{diag}\left(1-r_{-}^{\text{(log)}}(f)\right)A(f)Z,$		(23)
	$\displaystyle\nabla^{2}_{-}L_{\text{LSQ}}(f)$	$\displaystyle=Z^{\top}\operatorname{diag}\left(A(f)^{\top}r_{+}^{\text{(log)}}(f)\right)Z-Z^{\top}A(f)^{\top}\operatorname{diag}\left(r_{+}^{\text{(log)}}(f)\right)A(f)Z,$		(24)

where

$$r(f)=r_{+}(f)-r_{-}(f),\;\;r_{+}(f)=\max\left[r(f),0\right]\text{ and }r_{-}(f)=\max\left[-r(f),0\right],$$

and similarly

$$r^{\text{(log)}}(f)=r^{\text{(log)}}_{+}(f)-r^{\text{(log)}}_{-}(f),\;\;r_{+}^{\text{(log)}}(f)=\max\left[r^{\text{(log)}}(f),0\right]\text{ and }r^{\text{(log)}}_{-}(f)=\max\left[-r^{\text{(log)}}(f),0\right].$$

To summarize, the expression for the convex local upper bound to the quadratic approximation in Eq. (12) is

$$Q\left(c,f_{0};f\right)=L\left(f_{0}\right)+\left(f-f_{0}\right)^{\top}\nabla_{f}L\left(f_{0}\right)+\frac{1}{2}\left(f-f_{0}\right)^{\top}\nabla^{2}_{+}L\left(f_{0}\right)\left(f-f_{0}\right),$$

(25)

where $\nabla^{2}_{+}$ is used instead of $\nabla^{2}_{f}$ in the quadratic term. Here $Q$ depends parametrically on the counts data $c$, through the function $L$ (see Eq. (12)), and the expansion center $f_{0}$. The gradients of $L$ at $f_{0}$ are obtained from Eqs. (13) and (14), and the Hessian upper bounds are available from Eqs. (21) and (23). Note that the quadratic expression in Eq. (25) is not necessarily an upper bound of the data discrepancy functions, even locally, because we bound only the quadratic expansion. We employ the convex function $Q\left(c,f_{0};f\right)$ combined with convex constraints to generate descent steps for the generic non-convex optimization problem specified in Eq. (5).

The generic convex optimization addressed in Ref. [25] is

$$x^{\star}=\operatorname*{arg\,min}_{x}\left\{F(Kx)+G(x)\right\},$$

(26)

where $F$ and $G$ are convex, possibly non-smooth, functions and $K$ is a matrix multiplying the vector $x$. The ability to handle non-smooth convex functions is key for addressing the convex constraints of Eq. (5). In the primal-dual picture this minimization is embedded in a larger saddle point problem

$$\min_{x}\max_{y}\left\{y^{\top}Kx-F^{*}(y)+G(x)\right\},$$

(27)

using the Legendre transform or convex conjugation

$$F^{*}(y)=\max_{x}\left\{x^{\top}y-F(x)\right\},$$

(28)

and the fact that

$$F(x)=F^{**}(x)=\max_{y}\left\{y^{\top}x-F^{*}(y)\right\}$$

(29)

if $F$ is a convex function. The CP primal-dual algorithm of interest solves Eq. (27) by iterating on the following steps

	$\displaystyle y^{(n+1)}$	$\displaystyle=\operatorname*{arg\,min}_{y^{\prime}}\left\{F^{*}(y^{\prime})+\frac{1}{2\sigma}\|y^{(n)}+\sigma K\bar{x}^{(n)}-y^{\prime}\|^{2}_{2}\right\}$		(30)
	$\displaystyle x^{(n+1)}$	$\displaystyle=\operatorname*{arg\,min}_{x^{\prime}}\left\{G(x^{\prime})+\frac{1}{2\tau}\|x^{(n)}-\tau K^{\top}y^{(n+1)}-x^{\prime}\|^{2}_{2}\right\}$		(31)
	$\displaystyle\bar{x}^{(n+1)}$	$\displaystyle=2x^{(n+1)}-x^{(n)},$		(32)

where $n$ is the iteration index; $\sigma>0$ and $\tau>0$ are the primal and dual step sizes, respectively, and these step sizes must satisfy the inequality

$$\sigma\tau<\frac{1}{\|K\|^{2}_{2}}$$

where $\|K\|_{2}$ is the largest singular value of $K$. Because this algorithm solves the saddle point problem, Eq. (27), one obtains the solution to the primal problem, Eq. (26) along with its Fenchel dual

$$y^{\star}=\operatorname*{arg\,max}_{x}\left\{-F^{*}(y)-G^{*}(-K^{\top}y)\right\}.$$

(33)

The fact that both Eqs. (26) and (33) are solved simultaneously provides a convergence check: the primal-dual gap, the difference between the objective functions of Eqs. (26) and (33), tends to zero as the iteration number increases.

In some settings, the requirement $\sigma\tau<\frac{1}{\|K\|^{2}_{2}}$ may be impractical or too conservative, and the CP algorithm can instead be implemented with diagonal matrices $\Sigma$ and $T$ in place of $\sigma$ and $\tau$ [26], with the condition $\|\Sigma^{1/2}KT^{1/2}\|<1$ and the revised steps

	$\displaystyle y^{(n+1)}$	$\displaystyle=\operatorname*{arg\,min}_{y^{\prime}}\left\{F^{*}(y^{\prime})+\frac{1}{2}\|y^{(n)}+\Sigma K\bar{x}^{(n)}-y^{\prime}\|^{2}_{\Sigma^{-1}}\right\}$		(34)
	$\displaystyle x^{(n+1)}$	$\displaystyle=\operatorname*{arg\,min}_{x^{\prime}}\left\{G(x^{\prime})+\frac{1}{2}\|x^{(n)}-TK^{\top}y^{(n+1)}-x^{\prime}\|^{2}_{T^{-1}}\right\}$		(35)
	$\displaystyle\bar{x}^{(n+1)}$	$\displaystyle=2x^{(n+1)}-x^{(n)},$		(36)

where for a positive semidefinite matrix $A$ the norm $\|z\|_{A}$ is defined as $\sqrt{z^{\top}Az}$.

To apply the CP primal-dual algorithm to $Q$ for fixed $f_{0}$, we need to write Eq. (25) in the form of the objective function in Eq. (26). Manipulating the expression for $Q$ and dropping all terms that are constant with respect to $f$, we obtain

	$\displaystyle Q\left(c,f_{0};f\right)$	$\displaystyle=\frac{1}{2}f^{\top}K^{\top}(D-E)Kf-f^{\top}K^{\top}b,$		(37)
	$\displaystyle K$	$\displaystyle=\left(\begin{array}[]{c}K_{1}\\ K_{2}\end{array}\right)=\left(\begin{array}[]{c}A\left(f_{0}\right)Z\\ Z\end{array}\right),$		(42)
	$\displaystyle D$	$\displaystyle=\left(\begin{array}[]{cc}D_{1}&0\\ 0&D_{2}\end{array}\right),\;\;E=\left(\begin{array}[]{cc}E_{1}&0\\ 0&0\end{array}\right),\;\;b=\left(\begin{array}[]{c}b_{1}\\ b_{2}\end{array}\right)$		(49)
	$\displaystyle D_{1}$	$\displaystyle=\begin{cases}\operatorname{diag}\left[\hat{c}\left(f_{0}\right)\right]&\text{if }L=L_{\text{TPL}}\\ \mathbf{I}&\text{if }L=L_{\text{LSQ}}\end{cases},$
	$\displaystyle D_{2}$	$\displaystyle=\begin{cases}\operatorname{diag}\left[A\left(f_{0}\right)^{\top}r_{-}\left(f_{0}\right)\right]&\text{if }L=L_{\text{TPL}}\\ \operatorname{diag}\left[A\left(f_{0}\right)^{\top}r^{\text{(log)}}_{-}\left(f_{0}\right)\right]&\text{if }L=L_{\text{LSQ}}\end{cases},$
	$\displaystyle E_{1}$	$\displaystyle=\begin{cases}\operatorname{diag}\left[r_{-}\left(f_{0}\right)\right]&\text{if }L=L_{\text{TPL}}\\ \operatorname{diag}\left[r^{\text{(log)}}_{-}\left(f_{0}\right)\right]&\text{if }L=L_{\text{LSQ}}\end{cases},$
	$\displaystyle b_{1}$	$\displaystyle=\begin{cases}(D_{1}-E_{1})K_{1}f_{0}-r\left(f_{0}\right)&\text{if }L=L_{\text{TPL}}\\ (D_{1}-E_{1})K_{1}f_{0}-r^{\text{(log)}}\left(f_{0}\right)&\text{if }L=L_{\text{LSQ}}\end{cases},$
	$\displaystyle b_{2}$	$\displaystyle=D_{2}K_{2}f_{0}.$

The matrices $D$ and $E$ are nonnegative and depend on $c$ and $f_{0}$; $b$ is a vector which also depends on $c$ and $f_{0}$; and $K$ is a matrix that depends only on $f_{0}$. Both terms of $Q$ are functions of $Kf$ and accordingly $Q$ is identified with the function $F$ in the objective function of Eq. (26)

	$\displaystyle Q\left(c,f_{0};f\right)$	$\displaystyle=F_{Q}(Kf),$
	$\displaystyle F_{Q}(z)$	$\displaystyle=\frac{1}{2}z^{\top}(D-E)z-z^{\top}b.$		(50)

Because $Q$ is a convex function of $f$, $F_{Q}$ is convex as a function of $f$. The function $F_{Q}$, however, is not a convex function of $z$. Because $D$ and $E$ are non-negative matrices, $F_{Q}$ is a difference of convex functions of $z$,

	$\displaystyle F_{Q}(z)$	$\displaystyle=F_{Q+}(z)-F_{Q-}(z),$
	$\displaystyle F_{Q+}(z)$	$\displaystyle=\frac{1}{2}z^{\top}Dz-z^{\top}b,$
	$\displaystyle F_{Q-}(z)$	$\displaystyle=\frac{1}{2}z^{\top}Ez,$

where $F_{Q+}(z)$ and $F_{Q-}(z)$ are convex functions of $z$. That $F_{Q}(z)$ is not convex implies that $F_{Q}$ cannot be written as the convex conjugate of $F_{Q}^{*}(y)$, and performing the maximization over $y$ in Eq. (27) no longer yields Eq. (26).

To generalize the CP algorithm to allow the case of interest, we consider the function $F$ to be a convex-concave

$$F(z)=F_{+}(z)-F_{-}(z),$$

where $F_{+}$ and $F_{-}$ are both convex. The heuristic strategy for MOCCA is to employ a convex approximation to $F(z)$ in the neighborhood of a point $z=z_{0}$

$$F_{\text{convex}}(z_{0};z)=F_{+}(z)-z^{\top}\nabla F_{-}(z_{0});$$

(51)

(again we drop terms that are constant with respect to $z$). We then execute an iteration of the CP algorithm on the convex function $F_{\text{convex}}(z_{0};z)$; and then modify the point of convex approximation $z_{0}$ and repeat the iteration. The question then is how to choose $z_{0}$, the center for the convex approximation, in light of the fact that the optimization of $F$ in the CP algorithm happens in the dual space with $F^{*}$, see Eq. (30).

A corresponding primal point to a point in the dual space can be determined by selecting the maximizer of the objective function in the definition of the Legendre transform. Taking the gradient of the objective function in Eq. (29) and setting it to zero, yields

$$x=\nabla_{y}F^{*}(y).$$

(52)

We use this relation to find the expansion point for the primal objective function that mirrors the current value of the dual variables.

Incorporating the convex approximation $F_{\text{convex}}(z_{0};z)$ about the mirrored expansion point $z_{0}$ into the CP algorithm, yields the iteration steps for MOCCA

	$\displaystyle z^{(n+1)}_{0}$	$\displaystyle=\nabla_{y}F_{\text{convex}}^{*}\left(z^{(n)}_{0};y^{(n)}\right)=\Sigma^{-1}(y^{(n-1)}-y^{(n)}+\Sigma K\bar{f}^{(n-1)})$		(53)
	$\displaystyle y^{(n+1)}$	$\displaystyle=\operatorname*{arg\,min}_{y^{\prime}}\left\{F_{\text{convex}}^{*}\left(z^{(n+1)}_{0};y^{\prime}\right)+\frac{1}{2}\|y^{(n)}+\Sigma K\bar{f}^{(n)}-y^{\prime}\|^{2}_{\Sigma^{-1}}\right\}$		(54)
	$\displaystyle f^{(n+1)}$	$\displaystyle=\operatorname*{arg\,min}_{f^{\prime}}\left\{G(f^{\prime})+\frac{1}{2}\|f^{(n)}-TK^{\top}y^{(n+1)}-f^{\prime}\|^{2}_{T^{-1}}\right\}$		(55)
	$\displaystyle\bar{f}^{(n+1)}$	$\displaystyle=2f^{(n+1)}-f^{(n)},$		(56)

where $F_{\text{convex}}^{*}(z_{0};y)$ is convex conjugate to $F_{\text{convex}}(z_{0};z)$ with respect to the second argument; the first line obtains the mirror expansion point using Eq. (52) and the right hand side expression is found by setting to zero the gradient of the objective function in Eq. (54); the second line makes use of convex approximation $F_{\text{convex}}$ in the form of its convex conjugate; and the remaining two lines are the same as the those of the CP algorithm. For the simulations in this article, all variables are initialized to zero. Convergence of MOCCA, the algorithm specified by Eqs. (53) - (56), is investigated in an accompanying paper [30], which also develops the algorithm for a more general setting.

The MOCCA algorithm handles a fixed convex-concave function $F$, convex function $G$, and linear transform $K$. In order to apply it to the spectral CT data fidelity, we propose: employing the local quadratic expansion in Eq. (37) to which we apply MOCCA, re-expand the spectral CT data discrepancy at the current estimate of the material maps, and iterate this procedure until convergence. We refer to iterations of the core MOCCA algorithm as “inner” iterations, and the process of iteratively re-expanding the data discrepancy and applying MOCCA are the “outer” iterations. Because MOCCA allows for non-smooth terms, the convex constraints described in Sec. II-B can be incorporated and the inner iterations aim at solving the intermediate problem

$$f^{*}=\operatorname*{arg\,min}_{f}\left\{\sum_{w\ell}Q\left(c_{w\ell},f_{0};f\right)+\sum_{i}\delta(P_{i})\right\}.$$

(57)

For the remainder of this section, for brevity, we drop the constraints and write the update steps taking only for the spectral CT data fidelity. The full algorithm with the convex constraints discussed in Sec. II-B can be derived using the methods described in [27] and an algorithm instance with TV constraints on the material maps is covered in Appendix C.

In applying MOCCA to $Q\left(c_{w\ell},f_{0};f\right)$, we use the convex and concave components from $F_{Q}$ in Eq. (50) to form the local convex quadratic expansion needed in MOCCA, see Eq. (51),

$$F_{Q,\text{convex}}(z)=\frac{1}{2}z^{\top}Dz-(z-z_{0})^{\top}(b+Ez_{0}).$$

(58)

The corresponding dual function

$$F_{Q,\text{convex}}^{*}(y)=\frac{1}{2D}\|y+b+Ez_{0}\|_{2}^{2}+z^{\top}_{0}Ez_{0},$$

(59)

is needed to derive the MOCCA dual update step at Eq. (54). We note that because the material maps $f$ enter $Q\left(c_{w\ell},f_{0};f\right)$ only after linear transformation, $Kf$, and comparing with the generic optimization problem in Eq. (26), we have $G(f)=0$ for the present case where we only consider minimization of the data discrepancy.

In using an inner/outer iteration, a basic question is how accurately does the inner problem need to be solved. It turns out that it is sufficient to employ a single inner iteration, so that effectively the proposed algorithm no longer consists of nested iteration loops. Instead, the proposed algorithm performs re-expansion at every iteration:

	$\displaystyle f_{0}$	$\displaystyle=\bar{f}^{(n)}$		(60)
	$\displaystyle\Sigma^{(n)}$	$\displaystyle=|K_{1}(f_{0})|\mathbf{1}/\lambda;\;\;T^{(n)}=\lambda|K_{1}(f_{0})|^{\top}\mathbf{1}$		(61)
	$\displaystyle z^{(n+1)}_{0}$	$\displaystyle=(\Sigma^{(n)})^{-1}\left(y^{(n-1)}-y^{(n)}+\Sigma^{(n)}K_{1}(f_{0})\bar{f}^{(n-1)}\right)$		(62)
	$\displaystyle y^{(n+1)}$	$\displaystyle=(D_{1}(f_{0})+\Sigma^{(n)})^{-1}\left[D_{1}(f_{0})\left(y^{(n)}+\Sigma^{(n)}K_{1}(f_{0})\bar{f}^{(n)}\right)-\Sigma^{(n)}\left(b_{1}(f_{0})+E_{1}(f_{0})z^{(n+1)}_{0}\right)\right]$		(63)
	$\displaystyle f^{(n+1)}$	$\displaystyle=\bar{f}^{(n)}-T^{(n)}K_{1}(f_{0})^{\top}y^{(n+1)}$		(64)
	$\displaystyle\bar{f}^{(n+1)}$	$\displaystyle=2f^{(n+1)}-f^{(n)},$		(65)

where $f^{(0)}$, $f^{(-1)}$, $\bar{f}^{(0)}$, $y^{(0)}$, and $y^{(-1)}$ are initialized to zero vectors.

Before explaining each line of the one-step spectral CT algorithm specified by Eqs. (60)-(65), we point out important features of the use of re-expansion at every iteration: (1) There are no nested loops. (2) The size of the system of equations is significantly reduced; note that only the first matrix block of $K$, $D$, $E$, and $b$ (see Eq. (37) for their definition) appears in the steps of the algorithm. By re-expanding at every iteration the set of update steps for the second matrix block becomes trivial. (3) Re-expanding at every step is not guaranteed to converge, and an algorithm control parameter $\lambda$ is introduced that balances algorithm convergence rate against possible unstable iteration, see Sec. IV for a demonstration on how $\lambda$ impacts convergence. A similar strategy was used together with the CP algorithm in the use of non-convex image regularity norms, see [28].

The first line of the algorithm, Eq. (60), explicitly assigns the current material maps estimate to the new expansion point. In this way it is clear in the following steps whether $\bar{f}^{(n)}$ enters the equations through the re-expansion center or through the steps of MOCCA. For the spectral CT algorithm it is convenient to use the vector step-sizes $\Sigma^{(n)}$ and $T^{(n)}$, defined in Eq. (61), from the pre-conditioned form of the CP algorithm [26], because the linear transform $K_{1}(f_{0})$ is changing at each iteration as the expansion center changes. Computation of the vector step-sizes only involves single matrix-vector products of $|K_{1}(f_{0})|$ and $|K_{1}(f_{0})|^{\top}$ with a vector of ones, $\mathbf{1}$, as opposed to performing the power method on $K_{1}(f_{0})$ to find the scalar step-sizes $\sigma$ and $\tau$, which would render re-expansion at every iteration impractical. In Eq. (61), the parameter $\lambda$ enters in such a way that the product $\|\Sigma^{1/2}KT^{1/2}\|$ remains constant. For the preconditioned CP algorithm, $\lambda$ defined in this way will not violate the step-size condition. The dual and primal steps in Eqs. (63) and (64), respectively, are obtained by analytic computation of the minimizations in Eqs. (54) and (55) using Eq. (58) and $G(f)=0$, respectively. The primal step at Eq.(64) and the primal variable prediction step at Eq. (65) are identical to the corresponding CP algorithm steps at Eqs. (31) and (32), respectively. The presented algorithm accounts only for the spectral CT data fidelity optimization. For the full algorithm incorporating TV constraints used in the results section, see the pseudocode in Appendix C.