{'fraction_non_alphanumeric': 0.04067871032492786, 'fraction_numerical': 0.015870028854597918, 'mean_word_length': 4.353425117528543, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multi-Channel Convolutional Analysis Operator Learning for Dual-Energy CT Reconstruction

Objective. Dual-energy computed tomography (DECT) has the potential to improve contrast, reduce artifacts and the ability to perform material decomposition in advanced imaging applications. The increased number or measurements results with a higher radiation dose and it is therefore essential to reduce either number of projections per energy or the source X-ray intensity, but this makes tomographic reconstruction more ill-posed. Approach. We developed the multi-channel convolutional analysis operator learning (MCAOL) method to exploit common spatial features within attenuation images at different energies and we propose an optimization method which jointly reconstructs the attenuation images at low and high energies with a mixed norm regularization on the sparse features obtained by pre-trained convolutional filters through the convolutional analysis operator learning (CAOL) algorithm. Main results. Extensive experiments with simulated and real computed tomography (CT) data were performed to validate the effectiveness of the proposed methods and we reported increased reconstruction accuracy compared to CAOL and iterative methods with single and joint total-variation (TV) regularization. Significance. Qualitative and quantitative results on sparse-views and low-dose DECT demonstrate that the proposed MCAOL method outperforms both CAOL applied on each energy independently and several existing state-of-the-art model-based iterative reconstruction (MBIR) techniques, thus paving the way for dose reduction.


Introduction
Dual-energy computed tomography (DECT) is an energy-resolved X-ray imaging technique that offers two sets of attenuation measurements utilizing two different energy spectra acquired over the same anatomical region (Zhang, Liang, Wang, Cai, Wang, Tang,

Main Contribution
We extend the CAOL approach to multichannel settings and we develop a multi-channel convolutional analysis operator learning (MCAOL) framework that can exploit direct joint reconstruction, given the low-dose DECT measurements, where all the unknown images are reconstructed simultaneously by solving one combined optimization problem. This is the first time that MCAOL is applied to DECT image reconstruction and we demonstrate its superiority respect to CAOL. Furthermore, MCAOL requires considerably less memory compared to alternative DL approaches. The joint reconstruction approach is developed for a low-dose data acquisition protocol which consists of collecting data using a sparse angular sampling, using a different X-ray energy in consecutive steps and low X-ray photon counts.
In DECT, a reasonable prior assumption is that attenuation images at different energies can be expected to be structurally similar in the sense that an edge (e.g., an organ boundary) that is present at one energy, is likely to be at same location and alignment with the other energies as well, even though the contrast between materials will be different at each energy.
MCAOL technique reconstructs attenuation images from the projection data combined with multi-channel filters trained on a dataset of reconstructed images.
The central idea of MCAOL is to learn unsupervised DECT muti-channel convolutional dictionaries that can provide a joint sparse representation of the underlined images by jointly learning filters for the different energies: each atom not only carries individual information for each energy individually but also inter-energy information.
By reconstructing DECT images using MBIR techniques in conjunction with MCAOL, the multi-energy information can be optimally used by allowing the images to "talk to each other" during the reconstruction process through the learned joint dictionaries, reducing noise while preserving image resolution. In order to deal with the extreme low-dose scenario, we model the Poisson model exactly and we solve the image optimization problem by using approximated quasi-Newton method with constrained memory to achieve accurate joint reconstruction with limited computational complexity.

Notation and Paper Organization
The remainder of this paper is structured as follows: Section 2 briefly reviews the singlechannel CAOL algorithm as initially proposed in (Chun & Fessler 2019) while Section 3 presents the proposed approach to extend CAOL to multi-channel problems. In Section 4 we first describe the physical model of X-ray DECT from the continuous to discrete domain and the Poisson noise model; furthermore, we illustrate the proposed method to minimize the MCAOL reconstruction objective function using the exact low-dose DECT model. Finally, in Section 5 comprehensive results of MCAOL on a numerical extended cardiac-torso (XCAT) phantom (Segars, Mahesh, Beck, Frey & Tsui 2008) and clinical acquisition of a patient's full body are shown together with a comparison of its performance with state-of-the-art algorithms for both the sparse-views and low-dose DECT reconstruction.
Throughout the paper we adopt the following notations: matrices or discrete operators and column vectors are written respectively in capital and normal boldface type, i.e. A and a, to distinguish from scalars and continuous variables written in normal weight; [a] j denotes the j-th entry of a; '(·) ' denotes the transposition; ' ' denotes the 2-dimensional (2-D) convolution operator; given an image x ∈ R J and a filter d ∈ R P , P < J, d x represents the 2-D convolution obtained by zero padding and then truncating the boundary such that the output dimension is the same as the input dimension J; an image x ∈ R J is represented by a column vector for algebraic operations and by a √ J × √ J matrix for 2-D convolutions. Finally, we indicate the set of positive real numbers with the notation ]0, +∞[.

Learning Convolutional Regularizers for Image Reconstruction: CAOL
In this Section we review the foundation of CAOL for MBIR which is achieved by solving an optimization problem of the form where x ∈ R J is the 2-D image to reconstruct, y ∈ R I is the observed measurement, L is a data-fidelity term that incorporates the measurement model-generally taking the form of a NLL function-and R is a regularizer weighted by β > 0; I and J are respectfully the dimension of the measurement (number of detectors) and dimension of the image (number of pixels). The minimization is carried out with the help of iterative algorithms such as modified expectation-maximization (EM) for emission tomography (ET) (De Pierro 1995) or penalized weighted least-squares (PWLS) combined with separable paraboloidal surrogate (SPS) for CT (Elbakri & Fessler 2002). The regularizer R is designed such that the reconstructed imagex(y) has desired properties, such as smoothness and sparsity of the gradient. It can be also trained so that x(y) can be sparsely represented as a linear combination of basic elements, or atoms, regrouped in a dictionary.
We consider the CAOL approach (Chun & Fessler 2019) where the image is sparsely represented with convolutional kernels (filters). In the analysis model, the image is represented with "sparsifying" filters d k ∈ R P by the analysis operator where z k ∈ R J is a sparse feature image vector of the same dimension as the image x. The filters d k ∈ R P are vectorized images of dimension P J that are regrouped in a dictionary D = {d k } ∈ R P ×K .
Learning the dictionary D from a dataset of training images {x l ∈ R J : l = 1, . . . , L} corresponds to finding a collection of filters D = {d k } obtained by the following nonconvex optimization problem where C is the constraints on D = {d k } and with the training analysis objective function F a defined as where z l,k ∈ R J is the feature image associated to the training image x l and the filter d k , · 0 is the sparsity-promoting l 0 semi-norm defined for all z = [z 1 , . . . , z J ] ∈ R J as where 1 A : R → {0, 1} denotes the indicator function of a set A ⊂ R, which is defined as 1 A (ξ) = 1 if ξ ∈ A and 1 A (ξ) = 0 if ξ / ∈ A, and α > 0 is a weight balancing between accuracy and sparsity. In Chun & Fessler (2019) the filters are enforced to satisfy the tight-frame conditions, i.e., where I K is the K × K identity matrix, to promote filters diversity. The optimization problem (3) is solved by alternating between minimization in z and in D. The joint minimization in D, z can be achieved with a Block Proximal Extrapolated Gradient method using a Majorizer (BPEG-M) algorithm, which can be implemented using the CONVolutional Operator Learning Toolbox (CONVOLT) (Chun & Fessler 2019, Chun 2019. BPEG-M alternates between minimization in z and D. The minimization in z is obtained with a hard-thresholding operator T : for all a = [a 1 , . . . , a J ] ∈ R J and for all β > 0, which provides a global minimizer for z → 1 2 a − z 2 2 + β z 0 , in such a way that The minimization in D is achieved with a proximal mapping with orthogonality constraints (PMOC), which can be implemented with a singular value decomposition (Chun & Fessler 2019). Finally the regularizer R in the minimization problem (1) is derived from the learned filters D as 3. Multi-channel Convolutional Analysis Operator Learning MBIR can be generalized to multi-channel imaging. Assuming we wish to reconstruct two images x 1 , x 2 ∈ R J of the same "object" from two independent measurements y 1 ∈ R I 1 and y 2 ∈ R I 2 corresponding to two modalities, multi-channel MBIR can be achieved by using an iterative algorithm to solve min where L 1 and L 2 are the data-fidelity terms for x 1 and x 2 , R mc is a multi-channel regularizer and ρ 1 , ρ 2 > 0 are weights. R mc is designed to exploit the inference between the 2 channels x 1 and x 2 , for example to promote structural similarities as proposed in Ehrhardt et al. (2014).
MCAOL is achieved by solving the following optimization problem, given the training image set (x 1,l , x 2,l ) with the training cost function F mc defined as where γ 1 , γ 2 > 0 are weights and the semi-norm · 1,0 on R J × R J is defined for all z 1 = [z 1,1 , . . . , z 1,J ] ∈ R J and for all z 2 = [z 2,1 , . . . , z 2,J ] ∈ R J as · 1,0 denotes the l 1,0 norm. It promotes joint sparsity, i.e., with zero and non-zero values at the same locations, of image features in all the modalities, that are encoded by the multi-channel dictionary D 1 , D 2 .

X-ray CT Discrete Model
In this section, we describe the CT discrete physical measurement process with the spectrum of the X-ray source beams composed of two different energies. We consider the case of 2-D slice-by-slice imaging systems. For image reconstruction we assume that the continuous attenuation image µ e (r) which denotes the linear attenuation coefficient at position r ∈ R 2 and the energy level e = 1, 2, can be represented by a linear combination of basis functions {b j } associated to a discrete sampling on a √ J × √ J Cartesian grid, where x e,j > 0 for all j = 1, . . . , J and all e = 1, 2. The line integral becomes a summation: where ) dl is the contribution of the j-th pixel to the i-th ray. The system matrix A is constructed as an under-determined matrix of dimensions I × J where I = N d × N θ with N d and N θ being respectively the number of detectors and N θ and the number of angles (projections), and is defined as [A] i,j = a i,j , ∀ i = 1, . . . , I, ∀ j = 1, . . . , J. The spectral X-ray mathematical discrete model is based on the Beer's law which provides the X-ray intensity after transmission. The expected number of detected photonsȳ i,e is then redefined as a function of the discrete image x e asȳ where x e = [x e,1 , . . . , x e,J ] ∈ R J is the vector of attenuation coefficients at source energy e, S e is the mean photons flux at the e-th energy bin, as we assume a mono-energetic intensity, and η e,i ∈ R + is a known additive term representing the expected number of background events (primarily from scatter). In the case of normal exposure, the number of detected photons follows a Poisson distribution, i.e., and the measurements at each energy bin e = 1, 2 are stored in a vector y e = [y e,1 , . . . , y e,N d ·N θ ] .
Although monochromatic X-ray source does not usually hold for scanners in clinical practice, a common effective strategy consists of applying a polychromaticto-monochromatic source correction pre-processing step (Whiting, Massoumzadeh, Earl, O'Sullivan, Snyder & Williamson 2006), and in the rest of the paper we will therefore assume that we have a monoenergetic source or that it has already been appropriately corrected.

Low-Dose CT Reconstruction
In case of low X-ray dose, since the photons counts can be very limited, the Gaussian approximation is no longer applicable as the logarithm of the data cannot be computed. We therefore chose to perform sparse view CT reconstruction from the raw measurements (y 1 , y 2 ) by solving the minimization problem (10), with positivity constraints on (x 1 , x 2 ), using the Poisson NLL functions L 1 and L 2 defined as and the trained regularizer R mc derived from the learned filters (D 1 , D 2 ) as in (18). Therefore, substituting (23) and (18) into the minimization (10), we obtain the following explicit expression for the MCAOL DECT reconstruction problem: We solve the minimization problem (24) by the alternating estimation of the sparse feature images and the linear attenuation images {x e : e = 1, 2}. Given the current estimates of the sparse coefficients {z t k : k = 1, . . . , K}, the image update x t e at iteration t is obtained through the following minimization problem x t e = arg min with Φ e (x e ) = ρ e L e (x e , y e ) + γ e 2 K k=1 d e,k x e − z t e,k 2 2 .
In this work, we utilized a limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) algorithm (Nocedal & Wright 2006, Chapter 7) to solve (25). The L-BFGS iterative solver estimates x t e starting the previous iterate x t−1 e . We define the first estimate as x and where B t,(q) is an approximate inverse Hessian of Φ t e evaluated at x t,(q) e . The matrix/vector (26) is directly computed (without storing B t,(q) ) from the m previous iterates x t,(q−p) e , p = 0, . . . , m − 1. An approximate solution of the line-search

Algorithm 2: MCAOL Reconstruction Algorithm
Input: Initial images (x 0 1 , x 0 2 ), DECT learned filters D = (D 1 , D 2 ), joint sparsity weight γ = (γ 1 , γ 2 ), penalty weights ρ = (ρ 1 , ρ 2 ), DE sinogram y = (y 1 , y 2 ), system matrix A, intensities (S 1 , S 2 ) , #outer iterations N outer . Output: Reconstructed images (x 1 , x 2 ) for t = 0, . . . , N outer − 1 do Update sparse codes (in parallel) ; for k,. . . ,K do sub-problem is obtained by backtracking to match the Wolfe Conditions. The iterative scheme (i.e., w.r.t. q) is repeated until either a convergence criterion is met. A more detailed explanation can be found in Bousse, Courdurier,Émond, Thielemans, Hutton, Irarrazaval & Visvikis (2020). We utilized the implementation proposed in Zhu, Byrd, Lu & Nocedal (1997). We also used the L-BFGS algorithm to minimize L 1 (·, y 1 ) and L 2 (·, y 2 ) (without penalty) in order to obtain initial images x 0 1 and x 0 2 . The other part of the alternating scheme is to update the sparse features update z t e,k given the current estimate of x t e the sparse features update z t e,k . This step is achieved using the multi-channel thresholding operator defined in (16). The pseudo-code for MCAOL reconstruction algorithm is detailed in Algorithm 2.

Validation
We validated the proposed methods on two different DECT low-dose acquisition setup. In particular, we analyzed the case of sparse-view DECT reconstruction with normal photon dose and the case of extreme low-photon counts with increased number of views. By approximating the dose as the product of the number of views and photon counts, the latter case represents a more challenging scenario since the overall dose considered is lower than the sparse-view case. Our implementation was based on CONVOLT (Chun 2019).

Methods Used for Comparison
The objective of the simulations with sparse views and normal X-ray source intensity is to demonstrate that MCAOL achieves improved accuracy compared reconstructing each energy separately by solving (1) with the CAOL regularizer defined in (9) and with the edge-preserving TV regularizers, as well as simultaneously by solving (10) with the JTV regularizer, respectfully defined as and where N j denotes the 8 nearest neighboring pixels of pixel j and ω j,k are weights (ω j,k = 1 for axial neighbors and ω j,k = 1/ √ 2 for diagonal neighbors), and ε > 0 is a small real value to ensure differentiability. For each method, we used the L-BFGS solver to estimate x 1 and x 2 .
The experiment with extreme low-counts aims at demonstrating that considering a weighted least-squares approximation of the log-likelihood function no longer guarantees effective reconstruction results, instead the exact Poisson statistics should be accounted. This results in a degradation of the performance of CAOL when optimized through the PWLS solver while using the quasi-Newton solver L-BFGS leads to improved qualitative and quantitative results.

Methodology
All experiments were validated by generating the DECT measurements as in Eq. (21) and then running n = 20 Poisson noise instances as in Eq. (22) from a ground truth (GT) image x e = [x * e,1 , . . . , x * e,J ] ∈ R J , e = 1, 2. As performance metrics, we considered the mean absolute bias (AbsBias) error function defined as follows where x [n] e,j indicates the reconstructed linear attenuation coefficient at image pixel j from the n-th Poisson noise replicate, R is the spatial region of interest and N R is the number of pixels in the region R. Furthermore, we compute the standard deviation (STD) defined as e,j . In this work R corresponds to the nonnegative pixels region of x e and is the same for both energy levels.
The simulations were repeated for all the methods, for different values of the regularization parameters in the objective functions (1) and (10) in order to plot AbsBias / STD curves (Fig. 5, 9 and 11). The quality of the reconstruction is assessed by the proximity of the curve to the origin. Training and reconstruction were performed according to the below-described settings.
Training The optimization problem (13) is minimized using the BPEG-M algorithm (Chun & Fessler 2019) with normalized input dataset. To investigate the trade-off between accuracy and features sparsity, we tested (13) with different values of γ 1 = γ 2 with filter (d 1,k , d 2,k ) of dimension P = 49 and number of filters K = 49. For each simulation, we tuned γ 1 = γ 2 by testing different values to investigate the effect. For all the datasets, we have used a training of L = 25 images for each energy e. Regarding the BPEG-M algorithm, we set the tolerance value equal to 10 −4 and the maximum number of iterations to 3 × 10 3 . For the CAOL training algorithm, we have used the same settings as detailed for MCAOL except that we tuned a single regularization weight α in the optimization problem (4) for each separate energy channel.
Reconstruction MCAOL and JTV reconstructions were achieved by solving (10) with R mc defined as (18) and (28) respectively, while CAOL and TV reconstructions by solving (1) for each energy bin e = 1, 2 separately with R defined as (9) and (27) respectively. MCAOL and CAOL were achieved using N outer = 300 outer iterations while the inner image update is obtained using the L-BFGS algorithm with 300 iterations. The (γ 1 , γ 2 )-values and β-values were the same as for training. TV and JTV reconstructions were achieved with the L-BFGS algorithm with 300 iterations. The measurements were obtained from the GT images x e outside of the training set and the reconstructions were repeated for each noise instance n, for a range of (ρ 1 , ρ 2 )-values with ρ 1 = ρ 2 and for a range of β-values, in order to obtain AbsBias-versus-STD curves.
We performed sparse-views and low-dose experiments on a simulated XCAT phantom and clinical data to assess the potential of the method for medical practice as detailed below. The experiments were conducted with fixed X-ray dose amount, i.e., by selecting the number of angles and the X-ray source intensity, and we evaluated the quality of the linear attenuation images reconstructed with different methods, both qualitatively and quantitatively.

Results on XCAT Phantom
For the unsupervised MCAOL and CAOL training, the numerical data consists of 1-mm pixel-width 512 × 512 torso axial slice images generated from the XCAT phantom for 60 keV and 120 keV energies.
We utilized 20 slice pairs from the XCAT phantom, each pair consisting of a slice at E 1 = 60 KeV and a slice at E 2 = 120 KeV, to train the filters. An additional slice pair-not part of the training dataset-was used to generate the projection data (i.e, following (21) and (22)) for reconstruction, as detailed below. We used the MCAOL weights parameters γ 1 = γ 2 = 800 and the CAOL parameter α = 0.01.    Fig. 2(a)) and separate learning with CAOL ( Fig. 2(b)). From a qualitative point of view, we observe that MCAOL filter pairs (d 1 , d 2 ) look similar as the edges are identical in the 2 energy images, as opposed to the CAOL filters.
In order to generate the sparse-view DECT projection measurements (21), we modeled the projector A with a 2-mm full width at half maximum (FWHM) resolution parallel beam system and we used 1-mm pixel-width 406×406 GT torso axial-slice images with attenuation coefficients x 1 , x 2 at energies 120 keV (high) and 60 keV (low) which differs from the training examples. The simulation consisted on generating sparse-view sinograms with 406 detector pixels and 60 regularly spaced projection angles, where 360 • is the full view rotation. A monochromatic source withS e = 10 5 incident photons and 100 background events was used to generate each sinogram.
To support the statement that joint sparsity of allows both images to inform each other, which makes the estimation of z 1 and z 2 more robust, we show the estimated feature maps in Fig. 3 obtained using the XCAT data. By comparing the estimated  sparse feature maps z e,k for e = 1, 2 and k = 1, . . . , 49 for separate reconstructions (Fig. 3(a) and Fig. 3(b)) and joint reconstruction (Fig. 3(c) and Fig. 3(d)), there are no similarities between the feature maps obtained from separate reconstructions while the feature maps obtained from joint reconstruction have similar structures. Fig. 4 shows the XCAT GT and the reconstruction images for both 60 keV and 120 keV energies obtained by MCAOL and the other algorithms used for comparison. The images are obtained using the parameters which corresponds to the minimum AbsBias shown in Figs. 5(a) and 5(b). It is worth noting that MCAOL (Fig. 4(d)) manages to substantially reduce the noise as compared with CAOL (Fig. 4(c)). Fig. 5(a) and Fig. 5(b) show the AbsBias against the STD results respectively for low and high X-ray source energy. Among the methods used for comparison, TV promotes sparsity of the gradient, while JTV promotes joint sparsity of the 2 gradients and therefore are particularly well-suited for XCAT. Despite this observation, it is possible to show that the minimum AbsBias obtained by MCAOL outperforms all other algorithms, or in other words by fixing the STD, the AbsBias achieved by MCAOL is always lower while it is possible to claim that by fixing the AbsBias, the STD of MCAOL is reduced.

Results on Clinical Data
We utilized images reconstructed from data acquired on Philips IQon Spectral CT and reconstructed with a MBIR technique (Philips IQon Elite Spectral CT product specifications 2018). All patients provided signed permission for the use of their clinical data for scientific purposes and anonymous publication of data.
The experiment was conducted in a similar fashion as for the XCAT simulation. We selected 22 slice pairs from a full body patient scan with 0.902-mm pixel-width and 512 × 512 image size for the training dataset corresponding to thorax. The energies used in this study are 70 keV and 140 keV. An additional slice pair was used to generate the projection data for reconstruction, as detailed below.
The pair of trained filters (d 1,k , d 2,k ) obtained by both MCAOL and CAOL unsupervised learning is shown in Fig. 6; we used the parameters γ 1 = γ 2 = 10 4 and the CAOL parameter α = 10 −4 . To generate the sparse-view DECT measurements (21), we used the same geometrical and noise settings as for the XCAT simulation except that we used 451 detector pixels and 451 × 451 GT thorax images with attenuation coefficients x 1 , x 2 at energies 140 keV (high) and 70 keV (low) which differs from the training examples.
In Fig. 7 we show the estimated feature maps obtained using the clinical data. As already noted previously with the XCAT data, by comparing the estimated sparse feature maps z e,k for e = 1, 2 and k = 1, . . . , 49 for separate reconstructions (Fig. 7(a) and Fig. 7(b)) and joint reconstruction (Fig. 7(c) and Fig. 7(d)), there are no similarities between the feature maps obtained from separate reconstructions while the feature maps obtained from joint reconstruction have similar structures.  In Fig. 8 the GT image and the reconstruction images for both energies and the different methods are shown; it is worth noting that the MCAOL reconstruction ( Fig. 8(d)) is less noisy than the CAOL reconstruction (Fig. 8(c)).
Figs. 9(a) and 9(b) report the AbsBias versus the STD plots and we obtain a similar behavior compared to the XCAT simulations; although the relative distance of the AbsBias among the simulated algorithms is reduced, MCAOL is still outperforming all  other methods constantly either by fixing AbsBias or STD, while the performance of CAOL is improving and it is close to the JTV solution accuracy.

Results for Low-Dose DECT
We conducted DECT reconstruction on the same set of data as in Section 5.4 but with a different CT acquisition setup; we substantially decreased the initial photon counts tō S e = 10 3 (reduction of 2 orders of magnitude compared to the previous experiments) and we doubled the number of views to 120. By approximating the total delivered X-ray dose as the product of the photons intensity times the number of views, it turns out that this scenario is considerably more challenging in terms of ill-posed problem with a total dose reduction of 50 times. We used this simulation to prove that MCAOL returns a more accurate solution compared to other priors. Furthermore, we prove that despite the higher computational complexity to minimize the exact Poisson NLL in (24) compared to solving the problem with a weighted least-squares approximated NLL, i.e., PWLS data-fit cost function, MCAOL achieves substantial improved bias accuracy compared to the PWLS solution.
To perform these experiments, we used the same optimal learned convolutional filters as obtained by the MCAOL training procedure detailed in Section 5.4 and the GT images in Fig. 8(a). Fig. 10 show the reconstruction images for both energies and different methods; MCAOL (Fig. 10(c)) accurately reconstruct the image features compared to all other methods and it is confirmed that the PWLS model performs poorly (Fig. 10(f)).
Figs. 11(a) and 11(b) show either that MCAOL is consistently outperforming the other methods in terms of accuracy and variance and that the Poisson NLL formulation leads to a noticeable improvement compared to the PWLS formulation as it is indicated  by comparing CAOL and CAOL-PWLS.

Discussion
The bias-variance trade-off analysis of the estimation results over the regularization parameters confirms that MCAOL allows to achieve the minimum absolute bias compared to CAOL and other MBIR state-of-the-art methods and also reduce standard deviation. Furthermore, MCAOL has the benefit of requiring less memory respect to DL methods because of the convolutional structure of the trained filters.
The MCAOL algorithm can be potentially extended to multi-energies e = 1, . . . , E since the l 1,0 semi-norm in Eq. (15) can be defined for a set of E-vectorised feature maps z 1 , . . . , z E . Each z e , e = 1, . . . , E is a column vector of dimension J × 1. Then the joint l 1,0 semi-norm is defined as While the statistical noise tends to be higher in the multi-energy case, on each sub-band the contribute of the noise is reduced since the noise is split on more energy bands. Therefore, we believe that evaluating the joint norm, i.e., non-zeros elements in the feature vectors in overlapping positions for all energies at the same time, will reduce the degradation due to the increased overall noise. The MCAOL framework allows to utilize any mixed norms for the jointly sparse regularization and other norms, such as the l 2,1 -norm which as proposed by Degraux et al. (2017) for convolutional synthesis operator learning, may also be considered.
In our experiments we have considered the product between the X-ray source intensity and the number of projection angles as an empirical measure for the total transmitted X-ray dose. While this metric gives a good approximation of the dose, we consider the analysis of the standardized measure of radiation dose, i.e., CT dose index (CTDI), as well as the absorbed dose as a follow-up study.
We account the open problems of how to optimally select both the regularization norm and regularization parameter according to the dataset for future algorithm development.
Although this work focuses on the multi-channel imaging reconstruction problem, we believe that our proposed method can be utilized in conjunction to DECT to taskoriented material decomposition problems. In particular, while an approach would be to design a material decomposition module in the image space which takes as input the MCAOL reconstructed images, a more compelling strategy would be designing a direct approach from sinograms to material images through MCAOL.
Furthermore, MCAOL method can be exploited for other multimodal imaging application such as PET/CT and PET/MRI. In the multimodal case, given the different intensity range on each channel, a further analysis on how to choose the NLL weights γ 1 = γ 2 in (10) should be conducted to properly balancing the information coming from the different modalities.
Finally, from a learning point of view, MCAOL training can be seen as a multichannel single layer unsupervised convolutional autoencoder (Chun & Fessler 2019, Appendix A) which paves the way to extend this approach to deeper autoencoder architectures to capture more complex features such as textures. The analysis and comparison of the proposed MCAOL approach with other supervised deep learning approaches is planned as a follow-up study. It is important to stress that MCAOL inherits a precise mathematical derivation and therefore it should not be susceptible of instabilities in the reconstruction which have been proven to occur with deep learning methods (Antun, Renna, Poon, Adcock & Hansen 2020).
We consider these problems as future development of the proposed algorithm.

Conclusion
In this work, we have extended the convolutional analysis operator framework to multichannel imaging and we have applied and extensively analyzed the proposed method to the DECT application. The presented results show that by using the information coming from both energies and allowing the channels to "talk to each other" a more accurate solution of the reconstruction problem can be achieved together with a reduction of the noise in the estimate. The coupling between energies is encapsulated by using an l 1,0 sparse mixed norm in the MCAOL optimization problems both for training and reconstruction. We obtain consistently better performances across different DECT acquisition scenarios from sparse-views to low-dose photon counts.