Pengcheng Zhang

Shuhui Ren

Yi Liu

Zhiguo Gui

Hong Shangguan

Yanling Wang

Huazhong Shu

Yang Chen

A Total Variation Prior Unrolling Approach for Computed Tomography Reconstruction

Keywords: deep learning, computed tomography, unrolling approach, convolutional neural network, total variation, Chambolle-Pock algorithm 1

Background: With the rapid development of deep learning technology, deep neural networks can effectively enhance the performance of computed tomography (CT) reconstructions. One kind of commonly used method to construct CT reconstruction networks is to unroll the conventional iterative reconstruction (IR) methods to convolutional neural networks (CNNs). However, most unrolling methods primarily unroll the fidelity term of IR methods to CNNs, without unrolling the prior terms. The prior terms are always directly replaced by neural networks.

Purpose:

In conventional IR methods, the prior terms play a vital role in improving the visual quality of reconstructed images. Unrolling the hand-crafted prior terms to CNNs may provide a more specialized unrolling approach to further improve the performance of CT reconstruction. In this work, a primal-dual network (PD-Net) was proposed by unrolling both the data fidelity term and the total variation (TV) prior term, which effectively preserves the image edges and textures in the reconstructed images.

Methods: By further deriving the Chambolle-Pock (CP) algorithm instance for CT reconstruction, we discovered that the TV prior updates the reconstructed images with its divergences in each iteration of the solution process. Based on this discovery, CNNs were applied to yield the divergences of the feature maps for the reconstructed image generated in each iteration. Additionally, a loss function was applied to the predicted divergences of the reconstructed image to guarantee that the CNNs' results were the divergences of the corresponding feature maps in the iteration. In this manner, the proposed CNNs seem to play the same roles in the PD-Net as the TV prior in the iterative reconstruction methods. Thus, the TV prior in the CP algorithm instance can be directly unrolled to CNNs.

Results:

The datasets from the Low-Dose CT Image and Projection Data and the Piglet dataset were employed to assess the effectiveness of our proposed PD-Net. Compared with conventional CT reconstruction methods, our proposed method effectively preserves the structural and textural information in reference to ground truth.

Conclusions:

The experimental results show that our proposed PD-Net framework is feasible for the implementation of CT reconstruction tasks. Owing to the promising results yielded by our proposed neural network, this study is intended to inspire further development of unrolling approaches by enabling the direct unrolling of hand-crafted prior terms to CNNs.

a more tailor-made unrolling approach to improve the CT image quality.

Owing to the potential to perform edge preservation in reconstructed images, the TV prior has received considerable attentions in the field of CT reconstruction. [START_REF] Sidky | Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[END_REF][START_REF] Sidky | Convex optimization problem prototyping for image reconstruction in computed tomography with the Chambolle-Pock algorithm[END_REF][START_REF] Han | Algorithm-enabled low-dose micro-CT imaging[END_REF] In this work, to further improve IR methods based on the TV prior, we proposed a primal-dual network (PD-Net), which unrolls both the data fidelity and TV prior terms of the cost function to the CNNs. Owing to the convergence guarantee, the CP algorithm effectively solves the inverse problem based on the TV prior. [START_REF] Sidky | Convex optimization problem prototyping for image reconstruction in computed tomography with the Chambolle-Pock algorithm[END_REF] In this work, we provided a detailed derivation of the CP algorithm instance for CT reconstruction. After analyzing this derived instance, we discovered that the inclusion of the TV prior in the CP algorithm instance is to update the reconstructed images with its divergences in each iteration of the solution process. Based on this discovery, CNNs were applied to yield the divergences of the feature maps for the reconstructed image generated in each iteration. Additionally, a loss function was applied to the predicted divergences of the reconstruct image to ensure that the yielded results of the proposed CNNs are indeed the divergences of the corresponding feature maps in the iteration. In this manner, the proposed CNNs are founded to provide the same roles in the PD-Net as the TV prior in the IR methods. Thus, the TV prior in the CP algorithm instance can be directly unrolled to CNNs. Our PD-Net is constructed by unrolling both the data fidelity term and the TV prior in the CP algorithm instance. By doing so, the PD-Net effectively preserves the structural and textural information in the reconstruction results. The main contribution of our proposed PD-Net is a direct improvement to unrolling approaches by unrolling the TV prior in the IR methods to the CNNs.

METHODS

2.A. CP Algorithm Based on TV Prior

By incorporating a data-fidelity term and a penalty or prior term in the cost function, model-based iterative reconstruction methods can effectively restore the desired images from sinogram data. The TV prior plays a vital role in CT reconstruction owing to its robust performance in preserving image edges and textures. [START_REF] Sidky | Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[END_REF][START_REF] Sidky | Convex optimization problem prototyping for image reconstruction in computed tomography with the Chambolle-Pock algorithm[END_REF] Based on the TV prior, the cost function of the inverse problem for CT reconstruction is defined as:

2 2 1 1 min ( ) 2 u Au g u λ   - + ∂     , ( 1
)
where A is the system matrix, g are the sinogram data, and u is the reconstructed image obtained by solving this minimization problem. The first term is the data fidelity term, and the second is the TV prior. The regularization parameter λ controls the balance of these two terms. The TV prior is the 1  -norm of the gradient-magnitude image, which is calculated by the operator ∂ (in Appendix).

Many kinds of optimization algorithms [START_REF] He | Optimizing a parameterized plug-and-play ADMM for iterative low-dose CT reconstruction[END_REF] have been proposed to solve the minimization problem (1) and to obtain the desired images. The CP algorithm [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] effectively achieves a robust primal-dual scheme solution owing to its convergence guarantee. For a given primal minimization,

{ } min ( ) ( ) x F Kx G x + , (2)
the dual maximization is expressed as:

{ } * * max ( ) ( ) T y F y G K y - - - , (3)
ACCEPTED MANUSCRIPT / CLEAN COPY where x and y are the vectors in spaces X and Y , K is a linear transform from X to Y , the T superscript denotes matrix transposition, and F and G are two convex functions whose convex conjugations are * F and * G , respectively. The conjugate of a convex function can be computed using the Legendre transform. [START_REF] Sidky | Convex optimization problem prototyping for image reconstruction in computed tomography with the Chambolle-Pock algorithm[END_REF] In the solution process, the CP algorithm iteratively computes the primal minimization and dual maximization and achieves convergence when the duality gap between them is zero. [START_REF] Sidky | Convex optimization problem prototyping for image reconstruction in computed tomography with the Chambolle-Pock algorithm[END_REF] Given convex functions F and G , the primal minimization and the dual maximization can be analytically solved with the proximal mappings

* [ ] prox F σ and
[ ] prox G τ , as discussed in Ref. [4]. For the convex function H , the proximal mapping is calculated by the following minimization:

2 z z z prox H z H z σ σ ′   ′ -   ′ = +       . ( 4
)
The derivation of the proximal mappings As derived in Ref. [4] for the minimization problem (1) based on the TV prior, the convex functions F and G in the CP algorithm are:

2 1 1 ( , ) ( ) 2 F y z y g z λ = - + , (5)
( ) 0 G x = , (6)
where x u = , y Au = , and z u = ∂ . Terms ( )

1 ( ) /( 1 ) i i i p p Au g σ σ + ← + - + 6:
( )

1 ( )/m a x , i i i I i i q q u q u λ σ λ σ + ← + ∂ + ∂ 1 7: 1 1 1 ( ) T i i i i u u A p div q τ τ + + + ← - + 8: 1 1 1 ( ) i i i i u u u u θ + + + ← + - 9: 1 i i ← + 10: until i N ≥ ACCEPTED MANUSCRIPT / CLEAN COPY
I y g z prox F y z z σ σ λ σ λ   - =     +   1 , (7)
[ ]( ) prox G x x τ = , (8)
where I 1 is an image with all pixels set to 1. The corresponding CP pseudocode is shown in Alg. 1. As derived from Equation ( 7), the role played by the TV prior plays in the CP algorithm is to update vector i u with vector 1 ( )

i div q τ + .

2.B. PD-Net

The learned reconstruction scheme based on the CP algorithm was well-investigated for CT reconstruction in Ref. [13]. The Learned Primal-Dual reconstruction [START_REF] Adler | Learned primal-dual reconstruction[END_REF] replaces the proximal operators in the CP algorithm with neural networks whose parameters are learned from big data. However, this learned reconstruction model only focuses on the unfoldment of the data fidelity term in cost function (1) to CNNs without considering the prior term. It is well-known that good CT reconstruction performance can be achieved by incorporating an elaborately designed prior into the cost function of the inverse problem. [START_REF] Chen | LEARN: Learned experts' assessment-based reconstruction network for sparse-data CT[END_REF] In this work, to improve the performance of the learned reconstruction scheme, our PD-Net is constructed by unrolling both the data-fidelity term and the TV prior in minimization problem (1) to the CNNs.

As shown in Alg. 1, there are three proximal operators for the vectors i p (in the projection domain), i q (in the image domain), and i u (in the image domain), respectively. Because 0 =0 p , the formula for the vector 1 i p + in Line 5 of Alg. 1 can be rewritten as: ( ) ( )

( ) 2 , L A = ∂ ; 1/ L τ ← ; 1/ L σ ← ; 1 θ ← ; 0 i ← ;
1 1 0 1 i i k i k k p Au g σ σ + -+ = ← - +  6: ( ) 1 1 0 i k i i k i k m m k d u λ σ ξ -+ + = = ← ∂ ⋅  ∏ 7: 1 1 1
T i i i i u u A p d τ τ + + + ← - + 8: 1 1 1 ( ) i i i i u u u u θ + + + ← + - 9: 1 i i ← + 10: until i N ≥ ACCEPTED MANUSCRIPT / CLEAN COPY ( ) ( ) 1 1 0 1 i i k i k k p A u g σ σ + -+ = = - +  . (9)
Thus, the vector 1 i p + represents the sum of the difference between the forward projection of image i u and sinogram data g . As 0 =0 q , the proximal operator for the vector 1 i q + in Line 6 of Alg. 1 can be reformulated as:

1 1 0 i k i i k i k m m k q u λ σ ξ -+ + = = = ∂  ∏ , ( 10
)
where max( , )

i I i i q u ξ λ σ = + ∂ 1 .
Obviously, vector 1 i q + represents the sum of the gradients of images i u . The corresponding formula 1 ( ) i div q + in Alg. 1 can be redefined as:

( )

1 1 1 0 ( ) i k i i i k i k m m k d div q u λ σ ξ -+ + + = = = = ∂⋅  ∏ , ( 11
)
where the operator

T ∂⋅ = -∂ ∂ (in Appendix
) represents the calculation of the divergences for an image. Thus, 1 i d + is the sum of the weighted divergences of images i u . Given parameters i ξ , Alg. 1 can be rewritten as Alg.

2 without changing the final reconstruction results. As mentioned, the role played by TV prior in the CP algorithm instance is to update vector i u with vector 1 ( ) i div q + (i.e., vector ACCEPTED MANUSCRIPT / CLEAN COPY ( )

i i i u u u θ + + + -.
Thus, the inclusion of the TV prior in the cost function ( 1) is to update the vector 1 i u + with its weighted divergences 1 i d + in each iteration of the solution process. In this work, six kinds of CNNs were proposed to replace the proximal operators of the vectors i p , i d , and i u , unrolling both the data fidelity term and the TV prior in minimization problem (1) to the CNNs.

The proposed PD-Net is an unrolling approach. In this work, a fixed number 10 N = of iterations were unrolled to form the PD-Net for simplicity and better utility. For the deep-learning solution, multi-channel convolutions are employed to improve the feature extraction capabilities of neural networks. In each convolution layer, the effective information of a variable is usually represented by multi-channel feature maps.

However, for the iteration algorithm, only a single vector of each variable is passed to the next iteration, such as u + in Alg. 2. Thus, for conventional unrolling approaches, multi-channel feature maps must be reduced to a single feature map at the end of each iteration. This feature map reduction inevitably decreases the feature extraction capability of reconstruction networks. To address this drawback, we used a strategy similar to the Learned Primal-Dual [START_REF] Adler | Learned primal-dual reconstruction[END_REF] and passed all feature maps to the next iteration. In this work, the number of feature maps that persist between iterations were set to 5 The architecture of the proposed PD-Net is shown in Fig. 1 and the corresponding pseudocode is displayed in Alg. 3. The CNN p θ Γ (the blue block in Fig. 1(a)) is employed to replace the operator in Line 5 of Alg. A p + , feature maps i d are applied to restore feature maps

1 i p + , 1 i d + ,
1 i d + ′ in CNN d θ
Λ (the green block in Fig. 1(a)).

Because feature maps i d are the divergences of feature maps i u , the restored feature maps Θ . Thus, the differences

1 i d + of feature maps 1 i d +
′ and 1 i u + can be employed to represent the divergences of feature maps 1 i u + for the reconstructed image +1 i u in the current iteration. Moreover, a loss function is applied to the predicted divergences of the reconstructed image (Fig. 1), to ensure that feature maps

p θ , p θ , d θ , d θ , u
θ , and u θ , and their parameters contain parameters i ξ from Alg. 2. All these parameters are learned from big data. In this work, a ray-driven method and a pixel-driven method [START_REF] Zeng | Unmatched projector/backprojector pairs in an iterative reconstruction algorithm[END_REF] are employed to perform the forward projection i Au and the back-projection 1 T i A p + , respectively. There are two outputs 1 N u + and 1 N q + for the PD-Net. During training, the 2  -norm is utilized as the loss to minimize the differences between predicted images 1 N u + and actual images û , and the differences between the predicted divergences 1 N d + and the divergences d calculated from actual images û . Therefore, the total loss of the PD-Net can be defined as:

ACCEPTED MANUSCRIPT / CLEAN COPY 2 2 1 1 2 2 total N N L u u d d γ + + = - + - , ( 12
)
where γ is the weight parameter. In this work, this parameter is set to 0.5.

As displayed in Fig. 1, six kinds of CNNs employed are closely combined in the PD-Net, which is trained in an end-to-end manner. In this work, after analytically deriving the CP algorithm instance, we discovered that the TV prior contributes mainly to updating the reconstructed images by its divergences in each iteration of the solution process. Based on this discovery, two CNNs are proposed to yield the divergences of feature maps i u for the reconstructed image in each iteration. Additionally, to guarantee that the yielded results of these two CNNs are the same as the divergences of feature maps +1 i u during the iteration, a loss function is applied to the predicted divergences of the reconstructed image. Hence, the two proposed CNNs can play the same role in the PD-Net as the TV prior in the IR methods. Thus, the TV prior in the CP algorithm instance can be unrolled to CNNs. In this work, by unrolling the TV prior to the CNNs, our PD-Net effectively preserves the structural and textural information in the reconstruction results.

2.C. Datasets

The Low-Dose CT Image and Projection Data was employed to assess the effectiveness of our proposed PD-Net. The provided projection data were acquired from a helical trajectory but are not suitable for fan-beam CT reconstruction. In this work, we utilized a fan-beam geometry to simulated projection data in 2D CT scanning settings through the forward projection of the provided CT images. Chest CT images from 20 patients were selected, and we randomly chose 3,357 images from 18 patients for training and 100 images from the remaining two patients for testing. We performed the forward projections for the chosen images by equally sampling along the rotation angle direction. The linear detector has 768 bins with physical lengths of 1 mm. The distance is 1068.0 mm from the X-ray source to the detector arrays, and 595.0 mm from the X-ray source to the center of rotation. The datasets comprise simulated projection data and the corresponding chosen CT images.

These simulated projection data are employed to reconstruct the 512×512-pixel CT images. The physical height and width of each pixel are both 0.5859 mm. The simulated projection data are the noise-free data y , without considering the noise introduced during acquisition. To model a realistic clinical environment, Poisson noise and electronic noise [START_REF] Niu | Sparse-view X-ray CT reconstruction via total generalized variation regularization[END_REF] were added to the noise-free data y by ( )

0 2 0 ˆln exp( ) (0, ) I Poisson I Normal ε = -+ y y , ( 13
)
where 0 I is the number of X-ray photons employed to penetrate the object, and 2 ε denotes the variance used to simulate the electronic noise from the equipment measurement error. For the normal dose dataset, the projection data were generated under a 360-view scanning condition. The noise was added to the projection data with the parameters 6 0 10 I = and 2 10 ε = . [START_REF] Xia | MAGIC: manifold and graph integrative convolutional network for low-dose CT reconstruction[END_REF] Under the 360-view scanning condition, we also simulated the low-dose datasets of different dose levels, namely 15%, 10% and 5%. To add noise of the different levels to the datasets, parameter 0 I was set to

2.D. Implementation Details

The PD-Net model was implemented on the PyTorch deep learning toolbox. [START_REF] Paszke | Automatic differentiation in pytorch[END_REF] To learn the parameters of the PD-Net model from the dataset, the loss function (12) was optimized end-to-end using the Adam optimizer. [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] The initial learning rate was set to

(b1) (c1) (d1) (e1) (f1) (a2) (b2) (c2) (d2) (e2) (f2) (a3) (b3) (c3) (d3) (e3) (f3)
ACCEPTED MANUSCRIPT / CLEAN COPY Our proposed PD-Net was compared to several widely used reconstruction algorithms, including the FBP algorithm (Ram-Lak filter), TV algorithm, [START_REF] Sidky | Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[END_REF] DuDoTrans, [START_REF] Wang | DuDoTrans: Dual-domain transformer provides more attention for sinogram restoration in sparse-view CT reconstruction[END_REF] DIR-III, [START_REF] Su | Generalized deep iterative reconstruction for sparse-view CT imaging[END_REF] and Learned Primal-Dual. [START_REF] Adler | Learned primal-dual reconstruction[END_REF] The Tomographic Iterative GPU-based Reconstruction (TIGRE) toolbox 39 was employed to execute the TV 260 algorithm. The training settings of DuDoTrans, DIR-III, and Learned Primal-Dual were based on Refs. [37], [38] and [13], respectively. The normalized mean square error (NMSE), peak signal-to-noise ratio (PSNR), and ACCEPTED MANUSCRIPT / CLEAN COPY structural similarity index measure (SSIM) [START_REF] You | Structurally-sensitive multi-scale deep neural network for low-dose CT denoising[END_REF] were employed to quantitatively evaluate the performance of our proposed PD-Net and the comparison models.

) (b1) (c1) (d1) (e1) (a2) (b2) (c2) (d2) (e2) (a3) (b3) (c3) (d3) (e3) a1

RESULTS

265

3.A. Validation with Normal Dose Data

To demonstrate the effectiveness of our proposed PD-Net, we performed a qualitative comparison on the normal dose data simulated from the Low-Dose CT Image and Projection Data. Figure 2 presents the To further assess the ability of the proposed PD-Net, we performed a CT reconstruction using normal dose data simulated from the Piglet dataset. The reconstructed images of the different methods are shown in Fig. 5.

(b1) (c1) (d1) (e1) (f1) (a2) (b2) (c2) (d2) (e2) (f2) (a3) (b3) (c3) (d3) (e3) (f3) (a1) (b1) (c1) (d1) (e1) (a2) (b2) (c2) (d2) (e2) (a3) (b3) (c3) (d3) (e3)
Their performances are visually differentiated in Fig. 5 and over the three ROIs in Fig. 6. The corresponding residual images of the reconstruction results are displayed in Fig. 7. All reconstruction methods efficiently reconstruct tomography images from the projection data. Compared to the competing algorithms, PD-Net suppresses artifacts and maintains structural features better (Fig. 6), yielding the most similar reconstruction results to the ground truth (Fig. 7). According to human perceptual evaluations, our proposed PD-Net exhibits better visual quality with respect to authentic structural and textural information (Fig. 6). Consistent with the visual results, the proposed PD-Net achieved the highest PSNR and SSIM, and the lowest NMSE (Table II). It can be seen that our proposed PD-Net indisputably improved these three metrics.

3.B. Validation with Low Dose Data

To evaluate the performance of our proposed PD-Net for different dose levels, a dataset simulated with 10% III. It is clear that our proposed PD-Net obtains better numerical values than all competing methods.

3.C. Validation with Sparse Data

To further demonstrate the effectiveness of our proposed PD-Net, we also performed a CT reconstruction using sparse data. The reconstructed images obtained under the scanning condition of 120 views are demonstrated in Fig. 11, and a quantitative analysis of ROIs is provided in IV, the PD-Net scored the highest PSNR and SSIM, and the lowest NMSE.

3.D. Network Design Evaluation

In this section, we present an efficacy assessment of the combination chosen for the proposed PD-Net

(b1) (c1) (d1) (e1) (f1) (a2) (b2) (c2) (d2) (e2) (f2) (a3) (b3) (c3) (d3) (e3) (f3)
ACCEPTED MANUSCRIPT / CLEAN COPY supplementary file). We removed the BN and ReLU operators after each convolution and added a parametric ReLU (PReLU) after the first two convolutions. A feature map of the outputs of CNN were employed to construct the network architecture of the Primal-Dual-Fusion (Sup. 3). In these two kinds of CNNs, outputs 1 i p + and 1 i u + were obtained by extracting the effective information from the feature maps However, the Primal-Dual-Fusion also blurred some local textures. As shown in Fig. 14, the proposed PD-Net exhibited better visual quality and texture preservation than the other network designs. To quantitatively evaluate the statistical properties of images reconstructed by different network designs, we calculated the average scores and standard deviations of the PSNR, SSIM, and NMSE results for all testing phantom results.

) (b1) (c1) (d1) (e1) (a2) (b2) (c2) (d2) (e2) (a3) (b3) (c3) (d3) (e3) a1
(a1) (b1) (c1) (d1) (e1) (f1) (a2) (b2) (c2) (d2) (e2) (f2) (a3) (b3) (c3) (d3) (e3) (f3) (a1) (b1) (c1) (d1) (e1) (a2) (b2) (c2) (d2) (e2) (a3) (b3) (c3) (d3) (e3)
As shown in Table V, our proposed PD-Net scored the highest PSNR and SSIM, and the lowest NMSE. These promising reconstruction performances indicate the efficacy of the combination chosen for the proposed PD-Net framework.

3.E. Robustness Analysis

In clinical settings, scanning doses may vary based on the task at hand. Training multiple networks for different dose levels makes it clinically infeasible to apply most learning based algorithms. [START_REF] Antun | On instabilities of deep learning in image reconstruction-Does AI come at a cost?[END_REF][START_REF] Bhadra | On hallucinations in tomographic image reconstruction[END_REF][START_REF] Gilton | Model adaptation for inverse problems in imaging[END_REF] To verify the robustness of our proposed PD-Net, we utilized a model trained with 10% dose data to test the results at different dose levels: 15% and 5%. The chest reconstruction results are displayed in Fig. 15. The images in the first, third and fifth rows reflect the reconstructed images of the 5%, 10% and 15% dose cases, respectively. The images over ROIs are displayed in the second, fourth and sixth rows of ACCEPTED MANUSCRIPT / CLEAN COPY and overall quality using a five-point scale (1 = "Unacceptable", 2 = "Substandard", 3 = "Acceptable", 4 = "Above average", and 5 = "Excellent"). The subjective criterion of diagnostic acceptance is applied to evaluate whether the image quality of the reconstruction results are acceptable for clinical diagnosis. The quality scores 365 were reported as the average scores and standard deviations means. The quality scores of reconstructed images are shown in table VII. It can be seen that our proposed PD-Net achieved the highest scores. Here, the readers only examined the diagnozability in the reconstructed images and did not actually perform any detection/diagnosis on the reconstruction results.

3.F. Parameter Experiments 370

ACCEPTED MANUSCRIPT / CLEAN COPY

The number of feature maps for i p , i d , and i u that persist between iterations directly affect the visual quality of the reconstructed images. In this work, we present the quantitative measures of the reconstructions with respect to the different feature map numbers. For these reconstruction results, the boxplots of the three selected metrics (i.e., PNSR, SSIM, and NMSE) are displayed in Fig. 16. The proposed PD-Net method achieves the poorest image quality metrics, when all feature map numbers are set to 1. However, these metrics are improved with an increase in the number feature maps. The improvement is limited, when the number is larger than 5. Thus, the feature maps passed between iterations were set to 5 in this work.

The hyper-parameter λ controls the tradeoff between the two per-pixel 2  losses in Equation [START_REF] Wang | Deep learning for tomographic image reconstruction[END_REF]. To find the optimal value of this hyper-parameter, we performed parameter experiments using normal dose data. Fig. 17 displays the boxplots of PNSR, SSIM, and NMSE for the proposed PD-Net with weights of 0.1, 0.5, 1, 5 and 10. It can be seen that the proposed PD-Net method achieved the best image quality metrics, when λ equals 0.5. Thus, we used this setting in this work.

DISCUSSION

When including the TV prior in the cost function of an inverse problem, the IR methods yield promising reconstructed images with critical feature preservation that is useful for diagnosis. In this work, to further improve upon these methods, we developed a PD-Net unrolling approach in which both the data-fidelity term ACCEPTED MANUSCRIPT / CLEAN COPY and the TV prior of the cost function were unrolled to the CNNs. By deriving the CP algorithm instance for CT reconstruction, we discovered that the main contribution of including the TV prior in IR methods is to update the reconstructed images by its weighted divergences in each iteration of the solution process. Based on this discovery, CNNs were proposed to yield the divergences of feature maps for the reconstructed image in each iteration. Additionally, we applied a loss function to the predicted divergences, ensuring that the yielded results of the proposed CNNs were indeed the divergences of feature maps for the reconstructed image generated in the iteration. In this manner, the proposed CNNs play the same roles in the PD-Net as the TV prior in the IR methods. Thus, the TV prior in the CP algorithm instance can be unrolled directly to the CNNs in the PD-Net.

The experimental results (Figs. 2-13) show that our proposed PD-Net framework achieves robust and superior performances in texture preservation.

Compared to the Learned primal-dual model, [START_REF] Adler | Learned primal-dual reconstruction[END_REF] Conventional IR methods adopt an iteration technique to solve the inverse problem. In each iteration, the TV prior contributes to updating the reconstructed images with its weighted divergences. [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF] However, the iteration technique may blur image edges, especially those of inconspicuous edges. On one hand, the noise in the images may lead to incorrect estimations of image divergences in each iterate. On the other hand, the inappropriate selection of divergence weights may also reduce the visual quality of the reconstruction results. To overcome these drawbacks, a loss function was applied to the predicted divergences to avoid the incorrect estimation of the image divergences in the PD-Net. Moreover, divergence weights were learned from big data, reducing the blur of reconstructed images. For conventional unrolling approaches, IR methods directly unroll to neural networks for CT reconstruction.

In unrolling approaches, [START_REF] Adler | Learned primal-dual reconstruction[END_REF][START_REF] Chen | LEARN: Learned experts' assessment-based reconstruction network for sparse-data CT[END_REF][START_REF] Xia | MAGIC: manifold and graph integrative convolutional network for low-dose CT reconstruction[END_REF][START_REF] Liu | SGD-Net: efficient model-based deep learning with theoretical guarantees[END_REF][START_REF] Wang | ADMM-based deep reconstruction for limited-angle CT[END_REF][START_REF] He | Optimizing a parameterized plug-and-play ADMM for iterative low-dose CT reconstruction[END_REF][START_REF] Xiang | FISTA-Net: learning a fast iterative shrinkage thresholding network for inverse problems in imaging[END_REF] the number of feature maps that persist between the iterations is usually set to 1.

However, the reduction of persistent feature maps between the iterations may decrease the feature extraction capability of CT reconstruction networks. As observed in Fig. 10, the image quality metrics were the poorest when the feature map number was set to 1. Thus, to improve the visual quality of the reconstructed images, more feature maps were employed to transmit effective feature information between the iterations. As discussed in Ref. [45], the increase in the persistent feature maps may also improve the convergence rate of the reconstruction networks in the unrolling approaches.

The runtimes of the different reconstruction algorithms on normal dose data are shown in Table VIII. It can be seen that learning-based reconstruction algorithms take many hours for training. However, the testing time for a single image was dramatically reduced. When the model is trained, the three learning based approaches mentioned in this work will take far less execution time than is consumed by the TV algorithm, which facilitates clinical work flows. Owing to the calculations of CNN convolutions for image divergences, the computational time of our PD-Net was slightly longer than that of the Learned Primal-Dual algorithm. More effort is needed to reduce the computational time of the PD-Net to meet the requirements of clinical practice.

ACCEPTED MANUSCRIPT / CLEAN COPY

In this work, the projection data simulated from real CT images were employed for the experiments. Based on these simulated projection data, our proposed PD-Net yielded promising results with high image quality and quantitative metrics. However, the simulated projection data differed from the clinical CT data. In future research, we plan to further validate and assess the performance of our proposed PD-Net method on clinical data.

The robustness analysis showed that the proposed PD-Net is capable of reconstructing images from projection data with slightly different distributions. However, the PD-Net may not work correctly with the increase of the data distribution changes at test time. Future work will also focus on the robustness analysis for the larger data distribution changes.

CONCLUSIONS

In this work, we have developed a PD-Net method for the unrolling approaches applied to CT reconstructions. To enhance the feature extraction capability of the unrolling approaches, the TV prior of the IR methods was unrolled to the CNNs in the PD-Net. The experimental results demonstrate that our proposed PD-Net framework is feasible for the implementation of CT reconstruction tasks. Compared with conventional primal-dual networks, [START_REF] Adler | Learned primal-dual reconstruction[END_REF] our proposed PD-Net method effectively preserves the structural and textural information in reference to the ground truth. To further improve the unrolling approaches, future work will focus on the unfolding of the other sophisticated prior terms that are widely used in IR methods.

APPENDIX: THE FORMS OF THE SPATIAL GRADIENT AND THE SPATIAL DIVERGENCE

The explicit forms of the spatial gradient ∂ and the spatial divergence ∂ ⋅ in two dimensions are defined herein as follows. We use x I ∈ to refer to an M N × image and , i j

x to represent the ( , )

i j th pixel of x . To specify the linear transform ∂ , we introduce two gradient images c x I ∂ ∈ and r x I ∂ ∈ : The operator div is defined as T -∂ . With the Equation (A.3), the spatial divergence ∂ ⋅ becomes ( ) ( ) a) Author to whom correspondence should be addressed. Electronic mail: gzgtg@163.com (Z. Gui) and chenyang.list@seu.edu.cn (Y. Chen)

T x div x x ∂ ⋅ = ∂ = -∂ ∂ , ( ) , 1 , , , 1 ( ) ( )

  step for the CP algorithm to solve the minimization problem (1).

2

 2

λ 1 . 2 

 12 in Equation (5) are both convex functions. Thus, the function ( , ) F y z is convex because the sum of two convex functions is also convex. The Algorithm. Pseudocode for N steps of the 2 -TV CP algorithm. The constant L is the 2  -norm of the matrix K , σ and τ are non-negative CP algorithm parameters. The vectors

Algorithm. 2 . 2 

 22 New pseudocode for N steps of the 2 -TV CP algorithm.1:

FIG. 1 .

 1 FIG. 1. Network architecture of (a) our proposed PD-Net and (b) the CNN

  and

  1 i

  generated in the current iteration. The vectors +1 i u represent the feature maps of the reconstructed image +1 i u in the current iteration. Combined with vector 1 T i

Λ

  contain the more divergence information of the reconstructed image +1 i u than those of CNN u θ

  Line 7 of Alg. 2 and the operator in Line 8 are implemented by CNN u θ Ψ (the purple block in Fig. 1(a)), whose output 1 i u + is the reconstructed image of the current iteration. The outputs of the PD-Net are the reconstructed image 1 N u + and its divergence 1 N d + which is generated by CNN d θ Ω (Fig. 1(b)) with feature maps 1 N d + after the N th iteration. In this work, there are 3 convolutions per CNN. The pixel size and number of channels of each convolution are shown in Fig. 1. All convolutions are followed by batch normalization (BN) and a rectified linear unit (ReLU), apart from the last ones in CNNs p parameters of these six kinds of CNNs are

×

  , respectively.[START_REF] Xia | MAGIC: manifold and graph integrative convolutional network for low-dose CT reconstruction[END_REF] The sparse dataset was simulated under the 120-view scanning condition. The noise was added to the projection data with 6 0 10 I = . To further investigate the ability of the proposed PD-Net, we also performed experiments on the Piglet dataset, 34 which includes four different sub-datasets. We tested PD-Net on the sub-dataset of the normal dose CT, which contains 565 images for training and 190 images for testing. Using the same scanning conditions for the normal dose dataset simulated from the Low-Dose CT Image and Projection Data, we also generated a 250 normal dose dataset from the Piglet dataset.

4 10 -

 10 and was decreased by 20% after each epoch. The experiments were 255 performed on a workstation computer platform with an Intel® Core™ i7-9700K CPU, 32 GB RAM, and an NVIDIA RTX 3090 GPU.

FIG. 2 .

 2 FIG. 2. CT reconstruction results on normal dose data simulated from the Low-Dose CT Image and Projection Data. The display windows are [-1050, 1950] HU, [-250, 450] HU, and [-320, 480] HU for images in the three rows, respectively.

FIG. 4 .FIG. 5 .

 45 FIG. 4. Residual images of normal dose data reconstructions for the different methods in Fig. 2. The display window is [-10, 10] HU.TABLE I. Quantitative measures (average scores ± standard deviations) of the reconstructions on normal dose data simulated from the Low-Dose CT Image and Projection Data. Metrics TV DuDoTrans DIR-III Learned Primal-Dual PD-Net PSNR 32.5454±1.0149 34.9937±0.7886 35.1899±1.0018 34.5478±0.7318 37.7627±1.0817 SSIM 0.7801±0.0356 0.8575±0.0169 0.8693±0.0178 0.8230±0.0175 0.9071±0.0182 NMSE 0.0126±0.0015 0.0074±0.0018 0.0071±0.0020 0.0085±0.0027 0.0039±0.0008

FIG. 6 .FIG. 7 .

 67 FIG. 6. Zoomed parts from Fig. 5.

  ACCEPTED MANUSCRIPT / CLEAN COPY reconstructed images of the proposed PD-Net and competing methods. Three reconstruction cases from the testing dataset are shown in the three rows. For a better evaluation of the visual quality of the reconstructed images in Fig.2, the images over three regions of interest (ROIs) marked by the blue rectangles are zoomed-in and depicted in Fig.3. The corresponding residual images of the reconstruction results are shown in Fig.4. The TV algorithm's reconstructed images were shown to be contaminated by the piecewise constant effect (Figs.2(b) and 3(b)), and the images reconstructed by the DuDoTrans algorithm suffered from small streak artifacts (Figs.2(c) and 3(c)). As shown in Figs.2(d) and 3(d), the DIR-III algorithm works well with noise suppression and artifact removal but blurs some local structural features. The Learned Primal-Dual algorithm improved the visual quality of the reconstructed images (Figs. 2(e) and 3(e)). However, some tiny structures were lost in the reconstruction results. In Figs.2 and 3, the proposed PD-Net produced promising images in terms of texture preservation (Fig.3(f)). Upon visual inspection, our proposed PD-Net achieves better reconstruction accuracy than those of the competing algorithms (Fig.4).

  doses was utilized for training and testing. The reconstructed images are shown in Figs. 8-10. Because the radiation dose was severely reduced, the reconstruction results of the TV algorithm were seriously degraded with the loss of textural and structural information (Figs. 8(b) and 9(b)). In Figs. 8 and 9, the DuDoTrans, DIR-III and the Learned Primal-Dual algorithms efficiently suppressed the noise in the reconstructed images. However, the DuDoTrans results were contaminated by strong artifacts (Figs. 9(c1)-(c3)). The DIR-III algorithm blurred some crucial structural information (Figs. 9(d1)-(d3)). Additionally, the Learned Primal-Dual's results suffered from the loss of texture details (Figs. 9(e1)-(e3)). Our proposed algorithm exhibits better visual quality with respect to fine structural details, as shown in Figs. 9 and 10. For the low-dose data, the statistical quantitative results of the different reconstruction algorithms are shown in Table

Fig. 12 .

 12 The difference images of the reconstruction results are displayed in Fig. 13. As observed in Figs. 11(b) and 12(b), there are severe streak artifacts in the images reconstructed by the TV algorithm. The DuDoTrans and DIR-III algorithms greatly suppressed the streak artifacts in the reconstruction results (Figs. 11 and 12). However, there are still some tiny ACCEPTED MANUSCRIPT / CLEAN COPY streak artifacts in the DuDoTrans results (Figs. 11(c1)-(c3)). Figures 12(d1)-(d3) clearly show that the reconstruction results of the DIR-III algorithm suffered from losing diagnostic information. The Learned Primal-Dual algorithm effectively recovered the desired images from the projection data, but they suffered from the loss of structural information and spatial resolution (Figs. 12(e1)-(e3)). Compared with the Learned 315 Primal-Dual algorithm, our proposed PD-Net preserved more anatomical details, thereby yielding images of higher visual quality (Figs. 12(f1)-(f3)). Table IV lists the quantitative measures of the reconstructed images on sparse data. As shown in Table

FIG. 9 .

 9 FIG. 8. CT reconstruction results on low-dose data. The display windows are [-1050, 1950] HU, [-250, 450] HU, and [-320, 480] HU for images in the three rows, respectively.

FIG. 10 .FIG. 11 .

 1011 FIG. 10. Residual images of low-dose data reconstructions for the different methods in Fig. 8. The display window is [-10, 10] HU.TABLE III. Quantitative measures of the reconstructions on low-dose data. Metrics TV DuDoTrans DIR-III Learned Primal-Dual PD-Net PSNR 25.2906±1.4185 33.4349±1.0586 34.5303±0.9938 33.8402±1.0074 35.4082±1.1821 SSIM 0.5299±0.0484 0.8293±0.0227 0.8393±0.0249 0.7951±0.0285 0.8590±1.1821 NMSE 0.0673±0.1000 0.0105±0.0024 0.0082±0.0020 0.0097±0.0031 0.0067±0.0016

  respectively. Our proposed PD-Net further improved the Primal-Dual-Fusion by unrolling the 330 TV prior in the cost function via CNNs d images in the first row of Fig. 14 present the reconstruction results of the different network designs on normal dose data. The images over the ROIs of the reconstruction results are zoomed-in and depicted in the second row of Fig. 14. The Learned Primal-Dual suffered from the loss of subtle structural information (Figs.

FIG. 12 .

 12 FIG. 12. Zoomed parts from Fig. 11.

  ) and (b2)). The Primal-Dual-ReLU yielded reconstruction results with better subtle details than the Learned Primal-Dual (Figs. 14(c1) and (c2)). Compared with the Primal-Dual-ReLU, the Primal-Dual-Fusion further improved the reconstructed results by preserving more structural features (Figs. 14(d1) and (d2)).

Fig 15 .

 15 FIG. 14. CT reconstruction results on normal dose data for different network architectures. The display window is [-250, 450] HU.

FIG. 15 .

 15 FIG. 15. CT reconstruction results on normal dose data. The display window is [-250, 450] HU.

FIG. 16 .

 16 FIG. 16. Comparison of indices for the proposed method of different feature map numbers. Comparison of boxplots of (a) PSNR, (b) SSIM and (c) NMSE for reconstruction results.

  our proposed PD-Net improved the unrolling approaches by unrolling the TV prior in the IR methods to the CNNs. Moreover, the proposed CNNs in our PD-Net performed the same roles as the TV prior in the IR methods. The proposed PD-Net produced high-quality CT images with better anatomical feature preservations (Figs.2-13) and great improvements in image quality metrics is shown in Tables I-IV. The experimental results clearly demonstrate that unfolding the TV prior in the IR methods to the CNNs further improves the feature extraction capability of the unrolling approaches. ISTA-Net 28 is an unrolling approach, which is constructed by appropriately unrolling the iterative algorithm. To improve reconstruction performance, ISTA-Net adopts a CNN to replace the hand-crafted sparsity regularizer. Without giving the analytical formula of sparsity transform, ISTA-Net cannot fully unroll the prior term. Compared with ISTA-Net, the proposed PD-Net was constructed by unrolling both the data-fidelity term and the prior term. The PD-Net affords a tailor-made unrolling approach with the unfoldment of the hand-crafted TV prior.

2 )

 2 Using these two definitions, the spatial gradient ∂ can be written as ACCEPTED MANUSCRIPT / CLEAN COPY

TABLE II .

 II Quantitative measures of the reconstructions on normal dose data simulated from the Piglet dataset.

	Ground Truth	TV	DuDoTrans	DIR-III	Learned Primal-Dual	PD-Net
	(a1)
	Metrics	TV	DuDoTrans	DIR-III	Learned Primal-Dual	PD-Net
	PSNR 26.0673±1.3502 35.3099±1.2407 35.9514±1.2205	38.3275±1.2774	44.8953±1.6320
	SSIM	0.8736±0.0327	0.9041±0.0239	0.9131±0.0187		0.9276±0.0190	0.9952±0.0020
	NMSE	0.0285±0.0060	0.0057±0.0094	0.0030±0.0012		0.0019±0.0014	0.0004±0.0002

  Table I lists the average scores and standard deviations of PSNR, SSIM, and NMSE for all testing phantom results. Our proposed PD-Net achieved the best numerical values among all competing methods.

TABLE IV .

 IV Quantitative measures of the reconstructions on sparse data.

	Ground Truth	TV	DuDoTrans	DIR-III	Learned Primal-Dual	PD-Net
	Metrics	TV	DuDoTrans	DIR-III	Learned Primal-Dual	PD-Net
	PSNR 28.8094±1.4618 33.1754±0.9225 33.4723±0.7394	34.0128±0.8184	35.4554±1.0040
	SSIM	0.7241±0.0287	0.8114±0.0244	0.8201±0.0230	0.7937±0.0202	0.8498±0.0249
	NMSE	0.0316±0.0136	0.0113±0.0032	0.0104±0.0025	0.0092±0.0024	0.0066±0.0015

TABLE V .

 V Quantitative measures of the reconstructions for different network architectures.

	Ground truth	Learned Primal-Dual	Primal-Dual-ReLU	Primal-Dual-Fusion	PD-Net
	(a1)	(b1)	(c1)		(d1)	(e1)
		PSNR=36.0844	SSIM=0.8559 PSNR=37.1948	SSIM=0.9094 PSNR=38.3972	SSIM=0.9157 PSNR=38.6456	SSIM=0.9220
		(a2)	(b2)		(c2)	(d2)	(e2)
	Metrics Learned Primal-Dual Primal-Dual-ReLU Primal-Dual-Fusion	PD-Net
	PSNR	34.5478±0.7318	36.8346±1.0672	37.2678±1.0136	37.7627±1.0817
	SSIM	0.8230±0.0175	0.8976±0.0158	0.9010±0.0177	0.9071±0.0182
	NMSE	0.0085±0.0027	0.0048±0.0011	0.0041±0.0008	0.0039±0.0008

TABLE VI .

 VI Quantitative measures of the reconstruction results for robustness analysis.

	Dose Levels	Ground Truth	FBP		TV	PD-Net
	(a1)	(b1)			(c1)	(d1)
	5%	PSNR=14.6592	SSIM=0.2157 PSNR=23.6493	SSIM=0.4709 PSNR=34.7137	SSIM=0.8407
	(a2)	(b2) (b3)			(c2) (c3)	(d2) (d3)
	10%	PSNR=14.7070	SSIM=0.2333 PSNR=26.2779	SSIM=0.5445 PSNR=35.5817	SSIM=0.8565
		(b4)			(c4)	(d4)
		(b5)			(c5)	(d5)
	15%	PSNR=16.1496	SSIM=0.2620 PSNR=27.7702	SSIM=0.5928 PSNR=35.9654	SSIM=0.8641
		(b6)			(c6)	(d6)
		Metrics	5%		10%	15%
		PSNR 34.3973±1.2255 35.4082±1.1821 35.7857±1.1969
		SSIM	0.8425±0.0264	0.8590±1.1821	0.8662±0.0248
		NMSE	0.0084±0.0020	0.0067±0.0016	0.0061±0.0013

TABLE VII .

 VII Visual assessment scores by two radiologist readers.

	Algorithms	Ground truth	TV	DuDoTrans	DIR-III	Learned Primal-Dual	PD-Net
	Noise Reduction	-	2.90±0.8307	3.95±0.7703	4.00±0.7532	4.15±0.7263	4.55±0.5895
	Sharpness	-	3.25±0.6225	3.65±0.6427	3.70±0.6852	3.60±0.6633	4.05±0.5895
	Diagnostic Acceptability	-	3.55±0.4975	3.85±0.6825	3.90±0.6457	4.00±0.7746	4.50±05916
	Overall Quality	4.60±0.4899	3.15±0.5723	4.10±0.5308	4.15±0.4729	4.20±0.4000	4.55±0.4975

TABLE VIII .

 VIII Time costs of the different methods for training and testing on normal dose data.

		TV	DuDoTrans	DIR-III	Learned Primal-Dual	PD-Net
	Train	-	46 h	37 h	18 h	60 h
	Test	1963 s	0.2681 s	0.1647 s	0.1443 s	0.3715 s

A. Biguri, M. Dosanjh, S. Hancock, and M. Soleimani, "TIGRE: A MATLAB-GPU toolbox for CBCT image reconstruction," Biomed. Phys. Eng. Express 2, 055010 (2016).

ACKNOWLEDGMENTS

The authors would like to thank the institution, the Mayo Clinic, for providing the data (Low-Dose CT Image and Projection Data) used in this study. This work was supported in part by the State Key Project of Research and Development Plan under Grant 2022YFC2401600, in part by the National Natural Science Foundation of China under Grant T2225025, in part by the Key Research and development Programs in Jiangsu Province of China under Grants BE2021703 and BE2022768, in part by the Natural Science Foundation of Shanxi Province under Grant 201901D211246 and 202103021223308, and in part by the Research Project Supported by Shanxi Scholarship Council of China under Grant 2021-111.

CONFLICT OF INTEREST

The authors have no conflicts to disclose.


Shida Wang

Jalal Fadili

Peter Ochs

A Quasi-Newton Primal-Dual Algorithm with Line Search ⋆

Keywords: quasi-Newton, primal-dual algorithm, line search, saddlepoint problems, large scale optimization

Quasi-Newton methods refer to a class of algorithms at the interface between first and second order methods. They aim to progress as substantially as second order methods per iteration, while maintaining the computational complexity of first order methods. The approximation of second order information by first order derivatives can be expressed as adopting a variable metric, which for (limited memory) quasi-Newton methods is of type "identity ± low rank". This paper continues the effort to make these powerful methods available for non-smooth systems occurring, for example, in large scale Machine Learning applications by exploiting this special structure. We develop a line search variant of a recently introduced quasi-Newton primal-dual algorithm, which adds significant flexibility, admits larger steps per iteration, and circumvents the complicated precalculation of a certain operator norm. We prove convergence, including convergence rates, for our proposed method and outperform related algorithms in a large scale image deblurring application.

Introduction

In modern optimization, the datasets and dimensionality of the problems and parameters is vastly increasing. In the early stages of large scale optimization, the shift from second order to first order optimization could cope with the increasing dimensionality of the problems. Second order methods achieve a significant progress per iteration at the cost of a high computational load, since the computation of the second derivative (Hessian) and oftentimes its inverse are required, which is intractable in the large scale regime. In contrast, first order methods have a low computational effort per iteration but also less information about the objective. The gradient cannot capture curvature information and hence may fail to provide directions that allow for large steps. Nevertheless, for the large scale regime this exchange pays off.

However, the ever increasing dimensionality of the considered problems and datasets asks for faster algorithms. Motivated by classical optimization and the discussion above, algorithms at the interface of first and second order methods are key to reach the next level. Tractability in the (nowadays extremely) large scale regime requires methods that are as cheap as first order methods, while progressing as substantially as second order algorithms: Quasi-Newton methods. While they are known for their outstanding performance in unconstrained smooth optimization, their development for non-smooth (or constrained smooth) problems is insufficiently understood. As we discuss in related work below, most algorithmic development is either too simplistic, in the sense that only a diagonal metric is admitted, which can hardly capture second order information of the objective, or too theoretical, in the sense that a good performance is proved in theory while the implementation cost is on a par with that of second order methods. Algorithmic subproblems (e.g. the evaluation of the proximal mapping) that are easy to solve with respect to the Euclidean metric may become intractable with respect to another metric. We pursue the line of research initiated in [START_REF] Becker | A quasi-Newton proximal splitting method[END_REF][START_REF] Becker | On quasi-Newton forward-backward splitting: proximal calculus and convergence[END_REF] that considers both aspects as equally important. Key is the observation that quasi-Newton methods actually generate a metric of the specific type "identity ± low rank", which allows for an efficient proximal calculus (cf. Section 6) that unlocks the quasi-Newton power-well-known from classical optimization-in the area of optimization for Machine Learning. This idea was recently transferred to non-convex optimization in [START_REF] Kanzow | Globalized inexact proximal newton-type methods for nonconvex composite functions[END_REF][START_REF] Kanzow | Efficient regularized proximal quasi-Newton methods for large-scale nonconvex composite optimization problems[END_REF] and to monotone inclusion problems in [START_REF] Wang | Inertial quasi-newton methods for monotone inclusion: Efficient resolvent calculus and primal-dual methods[END_REF]. A special case of the latter setting comprises the extremely broad class of convex-concave saddle point problems, which has numerous applications in Machine Learning, Computer Vision, Image Processing and Statistics [START_REF] Chambolle | On the ergodic convergence rates of a first-order primal-dual algorithm[END_REF][START_REF] Chambolle | An introduction to continuous optimization for imaging[END_REF][START_REF] Goldstein | Adaptive primal-dual hybrid gradient methods for saddle-point problems[END_REF][START_REF] Valkonen | A primal-dual hybrid gradient method for nonlinear operators with applications to MRI[END_REF][START_REF] Applegate | Practical large-scale linear programming using primal-dual hybrid gradient[END_REF].

In this paper, we restrict our interest to saddle-point problems only. This focus allows us to design a quasi-Newton primal-dual algorithm that is tailored to this setting and therefore highly efficient and adaptable thanks to an additional line search procedure. This has several advantages as compared to a fixed step size: (i) the oftentimes expensive computation of the operator norm can be avoided, (ii) the choice of metric need not obey any static spectral restrictions, and (iii) in many situations larger steps and thus a faster convergence is observed.

Our main contribution is reduction of the gap between the outstanding performance of quasi-Newton methods in classical optimization and quasi-Newton methods for (non-smooth) convex-concave saddle point problems for modern optimization in Machine Learning. In detail, our contribution is the following:

1. We extend the line-search based primal-dual algorithm in [START_REF] Malitsky | A first-order primal-dual algorithm with linesearch[END_REF] (extension of [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] by line search) to incorporate a quasi-Newton metric with efficiently implementable proximal mapping; thereby aiming equally at theoretical convergence guarantees (including convergence rates) as well as highly efficient implementation.

2. We unlock the use of multi-memory quasi-Newton metrics (L-BFGS and SR1 method) via a compact representation for primal-dual algorithms, including their efficient implementation via a semi-smooth Newton solver. 3. The proposed algorithm outperforms the line-search based primal-dual algorithm (with identity metric) on a challenging image deblurring problem.

Related Work

Due to page limitations, for an extended discussion of quasi-Newton approaches in non-smooth optimization and the vast literature on primal-dual algorithms, we refer to [START_REF] Wang | Inertial quasi-newton methods for monotone inclusion: Efficient resolvent calculus and primal-dual methods[END_REF].

Non-smooth quasi-Newton. For a class of non-smooth problems that are given as a composition of a smooth function h and a non-smooth function g, [START_REF] Patrinos | Forward-backward truncated Newton methods for convex composite optimization[END_REF][START_REF] Stella | Forward-backward quasi-Newton methods for nonsmooth optimization problems[END_REF] combine quasi-Newton methods with forward-backward splitting via the forward-backward envelope. If g is an indicator function, [START_REF] Schmidt | Projected Newton-type methods in machine learning[END_REF][START_REF] Schmidt | Projected newton-type methods in machine learning[END_REF] proposed a projected quasi-Newton method which requires either solving a complicated subproblem or is restricted to a diagonal metric. Later, their work was extended by [START_REF] Lee | Proximal Newton-type methods for minimizing composite functions[END_REF] to a more general setting. [START_REF] Karimi | IMRO: A proximal quasi-Newton method for solving l1regularized least squares problems[END_REF][START_REF] Becker | A quasi-Newton proximal splitting method[END_REF][START_REF] Becker | On quasi-Newton forward-backward splitting: proximal calculus and convergence[END_REF] developed algorithms with efficient evaluation of the proximal operator with respect to a low-rank perturbed metric. It is worth to mention that in [START_REF] Becker | A quasi-Newton proximal splitting method[END_REF][START_REF] Becker | On quasi-Newton forward-backward splitting: proximal calculus and convergence[END_REF] the subproblem is a low dimensional root finding problem which can be solved efficiently. Inspired by [START_REF] Becker | On quasi-Newton forward-backward splitting: proximal calculus and convergence[END_REF], the work in [START_REF] Kanzow | Efficient regularized proximal quasi-Newton methods for large-scale nonconvex composite optimization problems[END_REF] applied the limited-memory quasi-Newton method on non-convex problems.

PDHG. Primal-Dual Hybrid Gradient (PDHG) is widely used to solve saddle point problems [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Chambolle | On the ergodic convergence rates of a first-order primal-dual algorithm[END_REF]. However, in order to guarantee the convergence of PDHG, the computation of the norm of a operator K is required. To avoid this disadvantage, [START_REF] Malitsky | A first-order primal-dual algorithm with linesearch[END_REF] combined line search with PDHG to get a new algorithm PDAL. For faster convergence, variable metric is being used [START_REF] Goldstein | Adaptive primal-dual hybrid gradient methods for saddle-point problems[END_REF]. It shows the potential of combining quasi-Newton methods and PDAL via a variable metric, however suffers again from the need to solve more complicated subproblems, which is remedied in [START_REF] Wang | Inertial quasi-newton methods for monotone inclusion: Efficient resolvent calculus and primal-dual methods[END_REF] for the more general class of monotone inclusion problems and hence builds the grounding for our proposed line search variant.

Problem Setup: A class of saddle point problems

Let X and Y be finite dimensional real vector spaces with inner product ⟨•, •⟩ and induced norm ∥ • ∥ = ⟨•, •⟩. We consider saddle point problems

min x∈X max y∈Y ⟨Kx, y⟩ + g(x) + h(x) -f * (y) , (1)
where g, h : X → R are proper, lower semi-continuous (lsc), convex functions with h, additionally, having an L-Lipschitz continuous gradient, f * : Y → R is a proper, lsc, convex function; the convex conjugate (Legendre-Fenchel conjugate) of a function f , and K : X → Y is a bounded linear operator with operator norm L K := ∥K∥ and adjoint K * . Moreover, throughout this paper we assume that

(1) has a saddle point. We remark that (1) is equivalent to the primal problem

min x∈X f (Kx) + g(x) + h(x) , (2)
and to the dual problem

max y∈Y -f * (y) + (g * ▽ h * )(-K * y) , (3)
where we use the fact that the conjugate of a sum of two function (g + h) * equals the infimal convolution of the conjugate functions g * ▽ h * .

Our quasi-Newton Primal-Dual Algorithm with Line Search

The primal-dual algorithm that we develop in this paper is an extension of [START_REF] Malitsky | A first-order primal-dual algorithm with linesearch[END_REF] to incorporate a variable metric of quasi-Newton type, which itself adds an efficient line search procedure to the primal-dual hybrid gradient (PDHG) algorithm [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] (aka Chambolle-Pock algorithm) and the extension in [START_REF] Lorenz | An inertial forward-backward algorithm for monotone inclusions[END_REF]. The handling of (a possibly) non-smooth functions essentially relies on evaluating the proximal mapping, which we define here for a function g and parameter τ with respect to an arbitrary metric M ∈ S α (X), α > 0. Here, S α is the set of bounded self-adjoint linear operators from Hilbert space X to X such that M -αI is positive semi-definite for each M ∈ S α . For simplicity, some notations are introduced: where I is the identity (Euclidean) metric. Algorithm 1 presents the proposed algorithm. The algorithm alternates between updates of the dual variable (4) and the primal variable [START_REF] Byrd | Representations of quasi-newton matrices and their use in limited memory methods[END_REF], where the line search is only implemented in the primal update. Let us discuss the algorithm for a fixed iteration k, i.e., we are given x k , y k-1 , σ k-1 , θ k-1 and a monotone decreasing sequence of β k . The discussion of the quasi-Newton type variable metric in Step (i) is deferred to Section 5. Step (ii) is a standard dual update step. In Step (iii), we perform the line search. We select a basic step size σk and, in the ith loop of the line search, we perform a trial step [START_REF] Byrd | Representations of quasi-newton matrices and their use in limited memory methods[END_REF] with the current σ k = σk • µ i and check if the breaking condition (6) is satisfied. If 'yes', the current iteration k is completed. If 'no', the new trial step size is reduced to σ k = σk • µ i+1 (in the subsequent (i + 1)th line search step). Here, M k ∈ S α (X) is symmetric positive definite and α ∈ (0, 1).If M k is chosen as an identity, then we recover the breaking (stopping) criterion used in [START_REF] Malitsky | A first-order primal-dual algorithm with linesearch[END_REF]. However, in this paper, we adopt a variable metric M k which is generated by quasi-Newton methods to exploit the local geometry of the function h. As a result, it is more likely to obtain a larger step size σ k and fewer inner loops for the line search procedure as compared to Algorithm 1 Quasi-Newton PDHG with Line Search Require:

∥x∥ 2 M = ⟨M x, x⟩ , x ∈ X .
[initial data] x 1 ∈ X, [initial data] y 0 ∈ Y , [maximal iteration count] N ≥ 0, [scaling of line search parameter] µ ∈ (0, 1), [extrapolation parameter] θ0, [initial dual step size] σ0, [tolerance weight] δ ∈ (0, 1), [primal-dual step ratio] +∞ > β ≥ β k+1 ≥ β k > 0, ∀k ∈ N.
Update for k = 1, . . . , N : (i) Compute M k according to a quasi-Newton framework (cf. Section 5).

(ii) Compute dual update step:

y k = prox σ k-1 f * (y k-1 + σ k-1 Kx k ) . ( 4
) (iii) Select σk ∈ [ β k-1 β k σ k-1 , (1 + θ k-1 ) β k-1 β k σ k-1
] and compute the quasi-Newton primal update step by: Line search: Find the smallest power i = 0, 1, 2, . . . such that

ȳk = y k + θ k (y k -y k-1 ) , x k+1 = prox M k τ k g x k -τ k M -1 k K * ȳk -τ k M -1 k ∇h(x k ) (5)
with

σ k = σk • µ i , θ k = σ k σ k-1 , and
τ k = β k σ k satisfy τ k σ k ∥Kx k+1 -Kx k ∥ 2 + 2τ k h(x k+1 ) -h(x k ) -∇h(x k ), x k+1 -x k ≤ δ∥x k+1 -x k ∥ 2 M k . (6)
End of for-loop the Euclidean version (M k = I). The employed metric is of type "identity ± low rank" for which the proximal mapping can be computed efficiently as shown in Section 6.

Remark 1. While the line search procedure is formulated for the primal problem, by duality, the primal problem can be interpreted as the dual of the dual problem and, thus, the dual problem as the primal problem. As a consequence, an equivalent algorithm with line search on the dual can be easily stated.

Discussion of computational cost for line search. In general, every loop of the line search procedure requires recomputing several quantities, including [START_REF] Byrd | Representations of quasi-newton matrices and their use in limited memory methods[END_REF] and Kx k+1 , h(x k+1 ) and ∥x k+1 -x k ∥ 2 M k in [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF]. While this seems to be expensive at first glance, often ( 6) is satisfied after 1-3 trial steps and hence large steps are taken with a low cost, as we underline in our experiments in Section 7. Nevertheless, the cost can be further reduced significantly in certain special cases, observed in [START_REF] Malitsky | A first-order primal-dual algorithm with linesearch[END_REF] and generalized here to our setting, whenever prox M k τ k g is a linear or affine operator.

1. If g(x) = ⟨c, x⟩, then prox M k τ k g (u) = u -τ k M -1 k c
and therefore, we obtain

x k+1 = prox M k τ k g (x k -τ k M -1 k K * ȳk -τ k M -1 k ∇h(x k )) = x k -τ k [M -1 k K * ȳk + M -1 k ∇h(x k ) + M -1 k c] . Kx k+1 = Kx k -τ k [KM -1 k K * ȳk + τ k KM -1 k ∇h(x k ) + τ k KM -1 k c] . 2. If g(x) = 1 2 ∥x -b∥ 2 , then prox M k τ k g (u) = (I + τ k M -1 k ) -1 [u + τ k M -1 k b]
and therefore, we obtain

x k+1 = prox M k τ k g (x k -τ k M -1 k K * ȳk -τ k M -1 k ∇h(x k )) = (I + τ k M -1 k ) -1 [x k -τ k M -1 k K * ȳk -τ k M -1 k ∇h(x k ) + τ k M -1 k b] , Kx k+1 = K(I + τ k M -1 k ) -1 [x k -τ k (M -1 k K * ȳk + M -1 k ∇h(x k ) -M -1 k b)] . 3. Let g(x) = δ H (x)
, where H refers to the hyperplane H := {u : ⟨u, a⟩ = b}.

Then

prox M k τ k g (u) = u + b-⟨u,a⟩ ∥a∥ 2 M -1 k M -1
k a and therefore, we obtain

x k+1 = prox M k τ k g (x k -τ k M -1 k K * ȳk -τ k M -1 k ∇h(x k )) = x k -τ k [M -1 k K * ȳk + M -1 k ∇h(x k )] + b -x k -τ k [M -1 k K * ȳk + M -1 k ∇h(x k )], a ∥a∥ 2 M -1 k M -1 k a , Kx k+1 = Kx k -τ k [KM -1 k K * ȳk + KM -1 k ∇h(x k )] + b -x k -τ k [M -1 k K * ȳk + τ k M -1 k ∇h(x k )], a ∥a∥ 2 M -1 k KM -1 k a .

Convergence Analysis of Algorithm 1

Let us now analyze the convergence of Algorithm 1. As for most variable metric primal-dual algorithms (cf. [START_REF] Davis | Convergence rate analysis of primal-dual splitting schemes[END_REF][START_REF] Combettes | Variable metric forward-backward splitting with applications to monotone inclusions in duality[END_REF][START_REF] Combettes | A forward-backward view of some primal-dual optimization methods in image recovery[END_REF][START_REF] Wang | Inertial quasi-newton methods for monotone inclusion: Efficient resolvent calculus and primal-dual methods[END_REF]), we require the following restriction for the change of the metric from one iteration to the next. Under this condition, we can generalize all convergence results from [START_REF] Malitsky | A first-order primal-dual algorithm with linesearch[END_REF] by adapting their proofs.

Assumption 1. Let α ∈ (0, +∞). (M k ) k∈N is a sequence in S α (X) such that ∃C M ∈ R, s.t.sup k∈N ||M k || ≤ C M < ∞ , (∃(η k ) k∈N ∈ ℓ 1 + (N))(∀k ∈ N) : (1 + η k )M k ⪰ M k+1 . (7)
Lemma 1. (i) There exists some σ > 0 such that σ k ≥ σ for any k ∈ N.

(ii) The line search terminates.

(iii) If β k ≡ β, θ k is bounded from above by some θ for any k ∈ N.

The proof is provided in Section B.1.

Theorem 1. Consider Problem (1) and let the sequence (x k , y k ) k∈N be generated by Algorithm 1 with β k ≡ β where Assumption 1 holds. Then (x k , y k ) k∈N is a bounded sequence and its cluster points are solutions of (1). Furthermore, if f * | domf * is continuous and σ k is bounded from above, then the whole sequence (x k , y k ) k∈N converges to a solution of (1).

The proof is provided in Section B.2.

We obtain the same ergodic convergence rate as in [START_REF] Malitsky | A first-order primal-dual algorithm with linesearch[END_REF], with respect to the primal-dual gap G x,ŷ which is the difference (gap) between the optimal primal objective value in (2) and the optimal dual objective value (3).

Theorem 2. Let the sequence (x k , y k ) k∈N be generated by Algorithm 1 with β k ≡ β where Assumption 1 holds and (x, ŷ) be some saddle point of (1). Then it holds for a constant

D = Π k∈N (1 + η k ) < +∞ that G x,ŷ ( XN , Ȳ N ) ≤ D s N 1 2β ∥x 1 -x∥ 2 M1 + 1 2 ∥y 1 -ŷ∥ 2 + σ 1 θ 1 D x,ŷ (y 0 ) = O 1 N , (8)
where

s N := N k=1 σ k , XN := N k=1 σ k x k+1 s N and Ȳ N := σ1θ1y 0 + N k=1 σ k ȳk σ1θ1+s N
. The last equality in (8) provides a simplified rate in terms of the big-O notation.

The proof is provided in Section B.3. Under the additional assumption that g is strongly convex, improved convergence rates can be derived when (β k ) k∈N is varied appropriately. Theorem 3. Assume g is γ-strongly convex and (x k , y k ) k∈N is generated by Algorithm 1 with

β k = β k-1 min{1 + γ C M β k-1 σ k-1 , C θ }
, ∀k ∈ N , and

β 0 > 0 , (9)
where C θ ∈ R + is a constant, Assumption 1 holds and (x, ŷ) be some saddle point of (1). Then, we have (θ k ) k∈N is bounded from above. Furthermore, we obtain

∥x N -x∥ = O(1/N ) and G x,ŷ ( XN , Ȳ N ) = O(1/N 2 ) ,
where ( XN , Ȳ N ) are the ergodic sequences defined in Theorem 2.

The proof is provided in Section B.4.

Remark 2. For the result in Theorem 3, δ = 1 is also admitted.

Computing and Representing the quasi-Newton Metric

In this section, we abuse notation in order to follow the conventions of quasi-Newton methods 3 . The metric M k is expected to be an approximation of the Hessian ∇ 2 h(x k ) at x k for the k-th iteration. The most popular quasi-Newton methods are the SR1 and BFGS methods (and their low-memory variants), which update M k by adding a rank-one modification (SR1 method)

M k+1 := M SR1 k+1 = M k + (y k -M k s k )(y k -M k s k ) ⊤ (y k -M k s k ) ⊤ s k , ( 10
)
or a rank-two modification (BFGS method)

M k+1 := M BF GS k+1 = M k + y k (y k ) ⊤ (s k ) ⊤ y k - M k s k (s k ) ⊤ M k (s k ) ⊤ M k s k , (11)
respectively, where

s k := x k+1 -x k and y k := ∇h(x k+1 ) -∇h(x k ) . (12)
In order to apply quasi-Newton methods on large-scale problems, m-limited memory quasi-Newton methods are adopted [START_REF] Kanzow | Efficient regularized proximal quasi-Newton methods for large-scale nonconvex composite optimization problems[END_REF], with the most popular version being L-BFGS [START_REF] Wright | Numerical optimization[END_REF]. It means that instead of generating M k via all previous s i and y i for i = 1, . . . , k and M 0 , for each k, the metric M k is re-computed based on M k,0 and the most recent m vectors s i and

y i for i = k -m + 1, . . . , k, if k ≥ m.
Usually, m is very small, such that only a small storage will be required. As pointed out by [START_REF] Byrd | Representations of quasi-newton matrices and their use in limited memory methods[END_REF], the matrices of the m-limited memory version of quasi-Newton methods have a compact representation of the form

M k = M k,0 + A k Q -1 k A ⊤ k , (13)
where M k,0 ∈ R n×n , n = dim(X), is a symmetric positive definite matrix, A k ∈ R n×m , and a symmetric and non-singular matrix

Q k ∈ R m×m (m ≪ n).
For limited memory BFGS (known as L-BFGS), we have the following blockmatrix representation k Y k , respectively. By using a spectral decomposition Q -1 = V ΛV ⊤ with orthogonal V ∈ R s×s and diagonal Λ ∈ R s×s , for some s ∈ N, (13) is transformed into the compact representation

A k = M k,0 S k Y k ∈ R n×2m and Q k = -S * k M k,0 S k -L k -L * k D k ∈ R 2m×2m , (14)
M k = M k,0 + U 1 U ⊤ 1 -U 2 U ⊤ 2 , (15)
for some U 1 ∈ R n×m and U 2 ∈ R n×m . In detail, since Λ is a diagonal matrix with eigenvalues λ i , i = 1, 2, . . . , s of Q -1 k on the diagonal, we decompose Λ = Λ 1 -Λ 2 where Λ 1 , given by (Λ 1 ) i,i = max(λ i , 0), i = 1, 2, . . . , s, corresponds to the positive eigenvalues and Λ 2 , given by (Λ 2 ) i,i = max(-λ i , 0), i = 1, 2, . . . , s, corresponds to the negative eigenvalues. In this way, we obtain

U 1 := (A k V )Λ 1/2 1 and U 2 := (A k V )Λ 1/2 2 . ( 16
)
Theoretically, it is guaranteed that

M k = M k,0 + U 1 U ⊤ 1 -U 2 U ⊤
2 is positive definite [START_REF] Fletcher | Practical methods of optimization[END_REF] if s k y k > 0 for any k ∈ N. However, in order to account for numerical rounding errors and the assumption that M k ∈ S α is bounded from above by some C M , we adopt a scaling version:

Mk = M k,0 + γ 1 U 1 U ⊤ 1 -γ 2 U 2 U ⊤ 2 , M k = min{ C M -α ∥ Mk ∥ 2 , 1} Mk + αI , (17)
where ∥ Mk ∥ 2 denotes the l 2 norm of matrix Mk and we set α = 0.01,γ 1 = 1, γ 2 = 1, C M = 50 in practice. There is an easy way to make sure that Assumption 1 is satisfied by setting

γ 1 = η k ∥U1∥ 2 and γ 2 = η k ∥U2∥ 2 with arbitrary η k ∈ ℓ 1 + .

Proximal Calculus and Efficient Implementation

The transformation in Section 5 enables us to compute the proximal mapping with respect to the metric in the form of ( 15) via the proximal calculus developed in [START_REF] Becker | On quasi-Newton forward-backward splitting: proximal calculus and convergence[END_REF], which we state here for completeness.

Theorem 4. Let B = B 0 + U 1 U ⊤ 1 -U 2 U ⊤ 2 ∈ S σ (R n ) with σ > 0, B 0 ∈ S σ (R n ) and U i ∈ R n×ri with rank r i (i = 1, 2). Set B 1 = B 0 +U 1 U ⊤ 1 .
Then, the following holds:

prox B g (x) = prox B0 g (x + B -1 1 U 2 α * 2 -B -1 0 U 1 α * 1 ) , (18)
where α * i , i = 1, 2, are the unique zeros of the coupled system L(α) = L(α 1 , α 2 ) = 0, where α = (α 1 , α 2 ) ∈ R r1+r2 and L = (L 1 , L 2 ) is defined by

L 1 (α 1 , α 2 ) = U ⊤ 1 (x + B -1 1 U 2 α 2 -prox B0 g (x + B -1 1 U 2 α 2 -B -1 0 U 1 α 1 )) + α 1 , L 2 (α 1 , α 2 ) = U ⊤ 2 (x -prox B0 g (x + B -1 1 U 2 α 2 -B -1 0 U 1 α 1 )) + α 2 . (19)
Here, L : R r1+r2 → R r1+r2 is Lipschitz continuous.

The computation of a possibly complicated proximal mapping prox B g (x) is reduced to a simple (by assumption) proximal mapping prox B0 g and a low dimensional root finding problem in [START_REF] Lorenz | An inertial forward-backward algorithm for monotone inclusions[END_REF], which we tackle by the semi-smooth Newton solver proposed in [START_REF] Becker | On quasi-Newton forward-backward splitting: proximal calculus and convergence[END_REF], formulated in Algorithm 2. It requires to solve Newtonlike equations where the classic Jacobian at the current iterate α k is replaced by the Clarke Jacobian ∂ c L(α k ) (see [START_REF] Clarke | Optimization and nonsmooth analysis[END_REF]), where we account for inexact solutions of these subproblems in terms of an error e k . For completeness, we also state the convergence result of [START_REF] Becker | On quasi-Newton forward-backward splitting: proximal calculus and convergence[END_REF] for Algorithm 2.

Theorem 5. If g is in addition a tame function, then the Lipschitz continuous function L is semi-smooth [START_REF] Bolte | Tame functions are semismooth[END_REF] and all elements of ∂ C L(α * ) are non-singular [START_REF] Becker | On quasi-Newton forward-backward splitting: proximal calculus and convergence[END_REF]. Therefore, if ρ k ≤ ρ, ∀k ∈ N, for some sufficiently small ρ and α 0 sufficiently close to α * , then the sequence generated by the Algorithm 2 is well-defined and converges to α * linearly. Additionally, if ρ k → 0, the convergence is superlinear.

Algorithm 2 Semi-smooth Newton method to solve L(α) = 0 in [START_REF] Lorenz | An inertial forward-backward algorithm for monotone inclusions[END_REF] Require:

[initial data] α0 ∈ R r , [maximum iterations] N Update for k = 0, • • • , N : (i) Select G k ∈ ∂ c L(α k ), compute α k+1 such that L(α k ) + G k (α k+1 -α k ) = e k ,
and e k ∈ R r is an error term satisfying

∥e k ∥ ≤ ρ k ∥G k ∥ and ρ k ≥ 0. (ii) if L(α k ) = 0 then terminate. End of for-loop
The tameness assumption is extremely mild, as it includes basically any function that occurs in practical applications, by excluding pathological special cases. For example this class of functions comprises all semi-algebraic functions [START_REF] Bolte | Tame functions are semismooth[END_REF].

Numerical Experiment

We apply our proposed algorithm for solving a challenging non-smooth image deblurring problem under a Poisson noise assumption [START_REF] Vardi | A statistical model for positron emission tomography[END_REF]. Given the observation b ∈ R nx×ny as n x × n y -sized image, the task is the following popular problem:

min x∈R nx ×ny + D KL (b, Ax) + γ∥Dx∥ 2,1 , (20)
which involves the Kullback-Leibler divergence as data fidelity measure h(x) := D KL (b, Ax) := i,j (Ax) i,j -b i,j log((Ax) i,j ) with respect to the blurred reconstruction Ax with known blur operator and a discrete total variation regularization term f (x) = γ∥Dx∥ 2,1 that is steered by a weight γ > 0 where D implements discrete spatial finite differences. We recast (20) into the saddle point problem:

min x∈R nx×ny max y∈R 2×nx ×ny ⟨Dx, y⟩ + δ R nx ×ny + (x) + D KL (b, Ax) -δ ∥•∥2,∞≤γ (y) (21)
and set g(x) := δ R nx ×ny + (x) and K = D in (1). Figure 1 compares several methods including PDHG with fixed stepsize (PDHG), PDHG with line search (PDAL), PDHG with fixed stepsize and variable metric (VarPDHG), PDHG with variable metric and line search (VarPDAL). The variable metric is generated by the limited memory BFGS method in Section 5. Figure 1 shows the primal gap where the optimal primal value was computed for 10000 iterations by running PDHG. For the update of the variable metric (17), we set γ 1 = 1.0 and γ 2 = 0.99. However, the Assumption 1 is not satisfied since it is not guaranteed by the construction of the metric that there is a sequence

(η k ) k ∈ ℓ 1 + (N) such that M k+1 ⪯ (1 + η k )M k .
Fortunately, we still observe convergence of PDHG with this variable metric. Figure 1 shows that a variable metric (VarPDHG) improves the convergence vs only using line search. However, our algorithm VarPDAL that combines both features is even faster, with the best performance when m = 9. Fig. 1: Performance evaluation for the experiment in [START_REF] Patrinos | Forward-backward truncated Newton methods for convex composite optimization[END_REF]. Algorithms are described in the text. Our algorithms, which combine a quasi-Newton variable metric with line search, outperform all other algorithms (that either use a variable metric VarPDHG or line search PDAL; or none of the two PDHG).

As a second experiment to test out accelerated version, we consider A = I and a constraint set C := {x|x ij ∈ [ϵ, 255]} since the grey value of each pixel should be less than 255 and be positive; We set ϵ = 0.1. In this case D KL (b, x) restricted on C is a strongly convex function. We apply the accelerated version of Algorithm 1 in Theorem 3, and we use the notation VarAPDAL. Similary, APDAL denotes the accelerated version of PDAL. The dashed line in 2 corresponds to O(1/N 2 ). We can observe that VarPDAL, APDAL can converge faster than O(1/N 2 ) as predicted by the convergence theorem.

In order to record the exact algorithmic details for our experiments, all code for experiments from this paper is available at https://github.com/ wsdxiaohao/VarPDAL.git. Fig. 2: Performance evaluation for the experiment for a strongly convex case. This figure reflects that our algorithms can retrieve O(1/N 2 ) convergence rate as PDAL, which is theoretically guaranteed by Theorem 3.

Conclusion

In this paper, we introduced a line search variant of a recently introduced quasi-Newton primal-dual algorithm. In contrast to related work, the employed quasi-Newton metric is of type "identity ± low rank", which captures significantly more second order information than a commonly used diagonal metric. We equally care for both, theoretical convergence guarantees including convergence rates as well as efficient practical implementation. The additional line search procedure usually leads to larger steps at a computational cost that pays off, which is confirmed by our numerical experiments.

A Preliminaries

There are several preliminaries we will use in the following section. The first one is a convergence result from [22, Lemma 2.2.2] of a special sequence which appears in B.2.

Lemma 2. Let a k ≥ 0 and let

a k+1 ≤ (1 + ν k )a k + ζ k , ν k ≥ 0, ζ k ≥ 0, k∈N ν k < ∞, k∈N ζ k < ∞. ( 22
)
Then, a k → A ≥ 0 for some A < +∞.

The following identity is called (cosine rule), which proves to be very useful.

2 ⟨a -b, c -a⟩ = ∥b -c∥ 2 -∥a -b∥ 2 -∥a -c∥ 2 ∀a, b, c ∈ X . ( 23
)
Another inequality appears many times in B.2 is the characteristic property of the proximal operator with respect to a symmetric positive definite matrix M :

x = prox M g (x) ⇐⇒ ⟨x -x, y -x⟩ M ≥ g(x) -g(y) ∀y ∈ X . (24)
If M = I is an identity matrix, then [START_REF] Schmidt | Projected newton-type methods in machine learning[END_REF] is the characteristic property of the standard proximal operator. Assume (x, ŷ) is a saddle point which solves (1).

Then we obtain

P x,ŷ (x) = g(x) + h(x) -g(x) -h(x) + ⟨K * ŷ, x -x⟩ ≥ 0 ∀x ∈ X , D x,ŷ (y) = f * (y) -f * (ŷ) -⟨K x, y -ŷ⟩ ≥ 0 ∀y ∈ Y , (25)
where P x,ŷ (x) and D x,ŷ (y) are convex. Then G x,ŷ (x, y) := P x,ŷ (x)+D x,ŷ (y) is the primal-dual gap. Without ambiguity, in the proofs, we may omit the subscript in P and D.

B Collection of Proofs

B.1 Proof of Lemma 1

It is a similar argument with the one in [START_REF] Malitsky | A first-order primal-dual algorithm with linesearch[END_REF].

(i)&(ii) σ k is decreased by µ ∈ (0, 1) and the inequality ( 6) is satisfied as long

as σ k < σ k := -1+ √ (4δα)/β k +1 2 L
where L = max{L, L K }. We introduce a

notation σ := -1+ √ (4δα)/β+1 2 L
. Since β k < β, we have σ k ≥ σ. We argument by induction. We assume σ 0 > µσ 0 and σ k-1 > µσ k-1 . For the case

σ k = σk , then σ k ≥ ( β k-1 β k )σ k-1 > µ( β k-1 β k )σ k-1 > µσ k > µσ. For the case σ k = µ i σk , σ ′ k = µ i-1 σk does not satisfy (6). It follows σ ′ k > σ k . Thus, σ k = µσ ′ k > µσ k ≥ µσ. (iii) By σ k ≤ σ k-1 1 + θ k-1 , we get θ k ≤ 1 + θ k-1 . Thus, θ k is bounded from above.

B.2 Proof of Theorem 1

The following proof is adapted from [START_REF] Malitsky | A first-order primal-dual algorithm with linesearch[END_REF]. Assume (x, ŷ) is a saddle point of problem 1 and β k ≡ β. By using [START_REF] Schmidt | Projected newton-type methods in machine learning[END_REF], we obtain the following two inequalities:

y k+1 -y k -σ k Kx k+1 , ŷ -y k+1 ≥ σ k (f * (y k+1 ) -f * (ŷ)) (26)
x k+1 -x k + τ k M -1 k K * ȳk + τ k M -1 k ∇h(x k ), x -x k+1 M k ≥ τ k (g(x k+1 ) -g(x)) (27) By using τ k = βσ k 1 β (x k+1 -x k ) + σ k M -1 k K * ȳk + σ k M -1 k ∇h(x k ), x -x k+1 M k ≥ σ k (g(x k+1 ) -g(x)) (28)
Similarly, we apply (24) on y k and obtain

y k -y k-1 -σ k-1 Kx k , y -y k ≥ σ k-1 (f * (y k ) -f * (y)) ∀y ∈ Y . (29)
Setting y = y k+1 and y = y k-1 respectively, we obtain

y k -y k-1 -σ k-1 Kx k , y k+1 -y k ≥ σ k-1 (f * (y k ) -f * (y k+1 )) ∀y ∈ Y , (30)
y k -y k-1 -σ k-1 Kx k , y k-1 -y k ≥ σ k-1 (f * (y k ) -f * (y k-1 )) ∀y ∈ Y . (31)
We deduce from (30) × θ k and θ k = σ k σ k-1 that:

θ k (y k -y k-1 ) -σ k Kx k , y k+1 -y k ≥ σ k (f * (y k ) -f * (y k+1 )) . (32)
By (31) × θ 2 k , we also get:

θ k (y k -y k-1 ) -σ k Kx k , θ k (y k-1 -y k ) ≥ σ k (θ k f * (y k ) -θ k f * (y k-1 )) . (33)
Summing (32) and (33) together, by using ȳk = y k + θ k (y k -y k-1 ), we obtain

ȳk -y k -σ k Kx k , y k+1 -ȳk ≥ σ k ((1 + θ k )f * (y k ) -θ k f * (y k-1 ) -f * (y k+1 )) .
(34) To sum up inequalties ( 26), ( 28) and (34) , we obtain

y k+1 -y k -σ k Kx k+1 , ŷ -y k+1 + 1 β (x k+1 -x k ) + σ k M -1 k K * ȳk + σ k M -1 k ∇h(x k ), x -x k+1 M k + ȳk -y k -σ k Kx k , y k+1 -ȳk ≥ σ k (f * (y k+1 ) -f * (ŷ)) + σ k (g(x k+1 ) -g(x)) + σ k ((1 + θ k )f * (y k ) -θ k f * (y k-1 ) -f * (y k+1 )) , (35)
Reorganizing the above inequality and using τ k = βσ k , we have

y k+1 -y k , ŷ -y k+1 + 1 β x k+1 -x k , x -x k+1 M k + ȳk -y k , y k+1 -ȳk + -σ k Kx k , y k+1 -ȳk + -σ k Kx k+1 , ŷ -y k+1 + σ k K * ȳk + σ k ∇h(x k ), x -x k+1 ≥ σ k (g(x k+1 ) -g(x)) + σ k ((1 + θ k )f * (y k ) -θ k f * (y k-1 ) -f * (ŷ)) , (36)
As in [START_REF] Malitsky | A first-order primal-dual algorithm with linesearch[END_REF], we still have:

-σ k Kx k , y k+1 -ȳk + -σ k Kx k+1 , ŷ -y k+1 + σ k K * ȳk , x -x k+1 = σ k Kx k -Kx k+1 , ȳk -y k+1 + σ k K x, ȳk -ŷ -σ k K * ŷ, x k+1 -x (37)
Adding σ k h(x k+1 ) -σ k h(x) on both sides of (36), we obtain:

y k+1 -y k , ŷ -y k+1 + 1 β x k+1 -x k , x -x k+1 M k + ȳk -y k , y k+1 -ȳk + -σ k Kx k , y k+1 -ȳk + -σ k Kx k+1 , ŷ -y k+1 + σ k K * ȳk + σ k ∇h(x k ), x -x k+1 + σ k h(x k+1 ) -σ k h(x) ≥ σ k (g(x k+1 ) -g(x) + (1 + θ k )f * (y k ) -θ k f * (y k-1 ) -f * (ŷ) + h(x k+1 ) -h(x)) . (38)
Combining (37) and (38), we have

y k+1 -y k , ŷ -y k+1 + 1 β x k+1 -x k , x -x k+1 M k + ȳk -y k , y k+1 -ȳk σ k Kx k -Kx k+1 , ȳk -y k+1 + σ k K x, ȳk -ŷ -σ k K * ŷ, x k+1 -x + σ k ∇h(x k ), x -x k+1 + σ k h(x k+1 ) -σ k h(x) ≥ σ k (g(x k+1 ) -g(x) + (1 + θ k )f * (y k ) -θ k f * (y k-1 ) -f * (ŷ) + h(x k+1 ) -h(x)) . (39)
By the definition of D(y) [START_REF] Stella | Forward-backward quasi-Newton methods for nonsmooth optimization problems[END_REF] and ȳk = y k + θ k (y k -y k-1 ), we have

(1 + θ k )f * (y k ) -θ k f * (y k-1 ) -f * (ŷ) -K x, ȳk -ŷ = (1 + θ k )(f * (y k ) -f * (ŷ) -K x, y k -ŷ ) -θ k (f * (y k-1 ) -f * (ŷ) -K x, y k-1 -ŷ ) = (1 + θ k )D(y k ) -θ k D(y k-1 ) . (40)
Using (40) and the definition of P (x), we deduce from (39) that

y k+1 -y k , ŷ -y k+1 + 1 β x k+1 -x k , x -x k+1 M k + ȳk -y k , y k+1 -ȳk + σ k Kx k -Kx k+1 , ȳk -y k+1 + σ k ∇h(x k ), x -x k+1 + σ k h(x k+1 ) -σ k h(x) ≥ σ k (P (x k+1 ) + (1 + θ k )D(y k ) -θ k D(y k-1 )) . (41)
From the line search condition (6), we have

σ k (h(x k+1 ) -h(x k ) -∇h(x k ), x k+1 -x k ) ≤ δ 2β ∥x k+1 -x k ∥ 2 M k - 1 2 σ 2 k ∥Kx k+1 -Kx k ∥ 2 . (42)
Additionally, by the convexity of h(x), we also have

h(x k ) -h(x) + ∇h(x k ), x -x k ≤ 0 . (43)
Combining (42) and σ k × (43), we get

σ k (h(x k+1 ) -h(x) -∇h(x k ), x k+1 -x ) ≤ δ 2β ∥x k+1 -x k ∥ 2 M k - 1 2 σ 2 k ∥Kx k+1 -Kx k ∥ 2 . (44)
Thus, it follows from (41) and (44) that

y k+1 -y k , ŷ -y k+1 + 1 β x k+1 -x k , x -x k+1 M k + ȳk -y k , y k+1 -ȳk + σ k Kx k -Kx k+1 , ȳk -y k+1 + δ 2β ∥x k+1 -x k ∥ 2 M k - 1 2 σ 2 k ∥Kx k+1 -Kx k ∥ 2 ≥ σ k (P (x k+1 ) + (1 + θ k )D(y k ) -θ k D(y k-1 )) . (45)
Using Cauchy-Schwarz inequality, we obtain

y k+1 -y k , ŷ -y k+1 + 1 β x k+1 -x k , x -x k+1 M k + ȳk -y k , y k+1 -ȳk + 1 2 σ 2 k ∥Kx k -Kx k+1 ∥ 2 + 1 2 ∥ȳ k -y k+1 ∥ 2 + δ 2β ∥x k+1 -x k ∥ 2 M k - 1 2 σ 2 k ∥Kx k+1 -Kx k ∥ 2 ≥ σ k (P (x k+1 ) + (1 + θ k )D(y k ) -θ k D(y k-1 )) . (46)
Applying [START_REF] Schmidt | Projected Newton-type methods in machine learning[END_REF], we deduce from (46) that

( 1 2 ∥y k -ŷ∥ 2 - 1 2 ∥y k+1 -y k ∥ 2 - 1 2 ∥ŷ -y k+1 ∥ 2 ) + ( 1 2β ∥x k -x∥ 2 M k - 1 2β ∥x k+1 -x k ∥ 2 M k - 1 2β ∥x -x k+1 ∥ 2 M k ) + ( 1 2 ∥y k -y k+1 ∥ 2 - 1 2 ∥ȳ k -y k ∥ 2 - 1 2 ∥y k+1 -ȳk ∥ 2 ) + 1 2 ∥ȳ k -y k+1 ∥ 2 + δ 2β ∥x k+1 -x k ∥ 2 M k ≥ σ k (P (x k+1 ) + (1 + θ k )D(y k ) -θ k D(y k-1 )) . (47)
Reorganizing the above inequalities, we obtain

1 2 ∥y k -ŷ∥ 2 + 1 2β ∥x k -x∥ 2 M k - 1 -δ 2β ∥x k+1 -x k ∥ 2 M k + σ k θ k D(y k-1 ) - 1 2 ∥ȳ k -y k ∥ 2 ≥ σ k (P (x k+1 ) + (1 + θ k )D(y k )) + 1 2 ∥ŷ -y k+1 ∥ 2 + 1 2β ∥x -x k+1 ∥ 2 M k . ( 48
) It follows from σk ≤ 1 + θ k-1 σ k-1 that σ k θ k ≤ σ 2 k σ k-1 ≤ σ2 k σ k-1 ≤ (1+θ k-1 )σ k-1 . Thus, 1 2 ∥y k -ŷ∥ 2 + 1 2β ∥x k -x∥ 2 M k - 1 -δ 2β ∥x k+1 -x k ∥ 2 M k + σ k-1 (1 + θ k-1 )D(y k-1 ) - 1 2 ∥ȳ k -y k ∥ 2 ≥ σ k (1 + θ k )D(y k ) + 1 2 ∥ŷ -y k+1 ∥ 2 + 1 2β ∥x -x k+1 ∥ 2 M k . (49)
Since (1 + η k )M k ⪰ M k+1 , we can obtain the following key inequality:

1 2 ∥y k -ŷ∥ 2 + 1 2β ∥x k -x∥ 2 M k - 1 -δ 2β ∥x k+1 -x k ∥ 2 M k + σ k-1 (1 + θ k-1 )D(y k-1 ) - 1 2 ∥ȳ k -y k ∥ 2 ≥ σ k (1 + θ k )D(y k ) + 1 2 ∥ŷ -y k+1 ∥ 2 + 1 2β(1 + η k ) ∥x -x k+1 ∥ 2 M k+1 . ( 50
) Set A k := 1 2 ∥y k -ŷ∥ 2 +σ k-1 (1+θ k-1 )D(y k-1 )+ 1 2β ∥x k -x∥ 2 M k .
Then, we deduce from (50) that

A k+1 ≤ (1 + η k )A k . ( 51
)
By Lemma 2, A k is bounded from above by some constant C. Thus, ∥y k -ŷ∥ and ∥x k -x∥ M k are both bounded. By the assumption that M k is uniformly bounded, ∥x k -x∥ is also bounded. As a result, we deduce from (50) that

k 1 -δ 2β ∥x k+1 -x k ∥ 2 M k + 1 2 ∥ȳ k -y k ∥ 2 ≤ k (1 + η k )A k -A k+1 ≤ C k η k + A 0 < +∞ . (52)
It implies that ∥x k+1 -x k ∥ M k → 0 and ∥ȳ k -y k ∥ → 0. So does ∥x k+1 -x k ∥ → 0, since (M k ) k∈N ⊂ S α (X). Since σ k > σ for some σ which is shown in Lemma 1 and β > 0 is fixed,

y k+1 -y k σ k = ȳk+1 -y k+1 σ k+1 → 0 as k → +∞ , ∥x k+1 -x k ∥ 2 M k τ k → 0 as k → +∞ . (53)
Since (x k , y k ) k∈N is bounded, we can extract a subsequence (x ki , y ki ) i∈N converging to some cluster point (x * , y * ). As in [START_REF] Malitsky | A first-order primal-dual algorithm with linesearch[END_REF], similarly, by using the lower semi-continuity of functions g and f * and the continuity of function h, we can pass the following two inequalities to the limit:

y ki+1 -y ki σ ki -Kx ki+1 , y -y ki+1 ≥ (f * (y ki+1 ) -f * (y)) ∀y ∈ Y , x ki+1 -x ki τ ki + M -1 ki K * ȳki + M -1 ki ∇h(x ki ), x -x ki+1 M k i = M ki (x ki+1 -x ki ) τ ki , x -x ki+1 + K * ȳki + ∇h(x ki ), x -x ki+1 ≥ (g(x ki+1 ) -g(x)) ∀x ∈ X . (54)
Thus, (x * , y * ) is the saddle point of (1). If, additionally, f * (y)| dom f * is continuous, then f * (y ki ) → f * (y * ) and D(y ki ) → 0 as i → +∞. From (50), we have 1 Π k j=1 (1+ηj ) A k is monotone. Setting x = x * and ŷ = y * in (50), by the boundedness of σ k and θ k , it follows that

lim k→∞ A k Π ∞ i=1 (1 + η i ) ≤ lim k→∞ A k Π k i=1 (1 + η i ) = lim i→∞ A ki Π ki j=1 (1 + η j ) ≤ lim i→∞ A ki = 0 (55) Since Π ∞ i=1 (1 + η i ) < +∞, we have lim k→+∞ A k → 0 which means x k → x * and y k → y * as k → +∞.

B.3 Proof of Theorem 2

We adapt the corresponding proof in [START_REF] Malitsky | A first-order primal-dual algorithm with linesearch[END_REF]. Let ϵ k := σ k P (x k+1 )+(1+θ k )D(y k )θ k D(y k-1 ) . Then we obtain the following inequality from (47),

1 2 ∥y k -ŷ∥ 2 - 1 2 ∥y k+1 -ŷ∥ 2 + 1 2β ∥x k -x∥ 2 M k - 1 2β ∥x k+1 -x∥ 2 M k - 1 2 ∥ȳ k -y k ∥ 2 ≥ ϵ k . (56)
By the assumption 1, we get

1 2 ∥y k -ŷ∥ 2 - 1 2 ∥y k+1 -ŷ∥ 2 + 1 2β ∥x k -x∥ 2 M k - 1 2β ∥x k+1 -x∥ 2 M k+1 (1 + η k ) - 1 2 ∥ȳ k -y k ∥ 2 ≥ ϵ k . (57) Since (1 + η k ) ≥ 1, it follows 1 2 ∥y k -ŷ∥ 2 - 1 2
∥y k+1 -ŷ∥ 2 (1 + η k ) + 1 2β ∥x k -x∥ 2 M k - 1 2β ∥x k+1 -x∥ 2 M k+1 (1 + η k ) ≥ ϵ k . ( 58
)
Let both sides of the above inequality be divided by Π k-1 i=1 (1 + η i ) and it is common to assume that an empty product yields identity i.e. Π 0 i=1 (1 + η i ) = 1. Thus,

1 2 ∥y k -ŷ∥ 2 Π k-1 i=1 (1 + η i ) - 1 2
∥y k+1 -ŷ∥ 2 Π k i=1 (1 + η i ) + 1 2β ∥x k -x∥ 2 M k Π k-1 i=1 (1 + η i ) - 1 2β ∥x k+1 -x∥ 2 M k+1 Π k i=1 (1 + η i ) ≥ ϵ k Π k i=1 (1 + η i ) . (59)
Summing up (59

) for k = 1, • • • , N , we obtain 1 2 ∥y 1 -ŷ∥ 2 + 1 2β ∥x 1 -x∥ 2 M1 ≥ N k=1 ϵ k Π k i=1 (1 + η i ) ≥ N k=1 ϵ k C . (60)
Here, we used the C = k∈N (1 + η k ) < +∞.

The following steps are similar with the ones in [START_REF] Malitsky | A first-order primal-dual algorithm with linesearch[END_REF].

N k=1 ϵ k =σ N (1 + θ N )D(y k ) + N k=2 [(1 + θ k-1 )σ k-1 -θ k σ k ]D(y k-1 ) -θ 1 σ 1 D(y 0 ) + N k=1 σ k P (x k+1 ) . (61)
Since D is convex,

σ N (1 + θ N )D(y N )+ N k=2 [(1 + θ k-1 )σ k-1 -θ k σ k ]D(y k-1 ) ≥ (σ 1 θ 1 + s N )D( σ 1 (1 + θ 1 )y 1 + N k=2 σ k ȳk σ 1 θ 1 + s N ) = (σ k θ 1 + s N )D( σ 1 θ 1 y 0 + N k=1 σ k ȳk σ 1 θ 1 + s N ) ≥ s N D( Ȳ N ) , (62)
where

s N = N k=1 σ k . Similarly, N k=1 σ k P (x k+1 ) ≥ s N P ( N k=1 σ k x k+1 s N ) = s N P ( XN ) . (63)
As a result,

G( XN , Ȳ N ) = P ( XN )+D( Ȳ N ) ≤ C s N 1 2β ∥x 1 -x∥ 2 M1 + 1 2 ∥y 1 -ŷ∥ 2 +σ 1 θ 1 D(y 0 ) . (64)

B.4 Proof of Theorem 3

The proof is also adapted from [START_REF] Malitsky | A first-order primal-dual algorithm with linesearch[END_REF]. From the update formula of β k , it follows that β k is decreasing. First, we are going to prove that θ k is bounded from above. It is not difficult but tedious. We know that if there exists a

C ∈ R + s.t θ k ≤ C 1 + θ k-1 then θ k is bounded. From this, it is sufficient to prove that β k-1
β k is uniformly bounded from above by some C θ . According to

β k = β k-1 min{1 + γ C M β k-1 σ k-1 , C θ } , ∀k ∈ N , and β 0 > 0 , (65)
we have that β k-1

β k = min{1 + γ C M β k-1 σ k-1 , C θ } ≤ C θ .
Second part, we are going to show the convergence rate. Since g is strongly convex, we obtain:

x k+1 -x k τ k + M -1 k K * ȳk + M -1 k ∇h(x k ), x -x k+1 M k ≥ (g(x k+1 ) -g(x)) + γ 2 ∥x k+1 -x∥ 2 . (66)
From Assumption 1, it follows that for any k ∈ N,

γ 2 ∥x k+1 -x∥ 2 ≥ γ 2C M ∥x k+1 -x∥ 2 M k+1 . (67)
Following the same way in which we got equation (48), by equation (66) and the assumption that

(1 + η k )M k ⪰ M k+1 , we obtain 1 2 ∥y k -ŷ∥ 2 - 1 2 ∥y k+1 -ŷ∥ 2 + 1 2β k ∥x k -x∥ 2 M k - 1 -δ 2β k ∥x k+1 -x k ∥ 2 M k - 1 2β k ∥x k+1 -x∥ 2 M k+1 (1 + η k ) - 1 2 ∥ȳ k -y k ∥ 2 ≥ ϵ k + γσ k 2 ∥x k+1 -x∥ 2 .
(68)

In order to obtain the following inequality, it is sufficient to assume δ ≤ 1. Thus,

1 2 ∥y k -ŷ∥ 2 - 1 2 ∥y k+1 -ŷ∥ 2 + 1 2β k ∥x k -x∥ 2 M k - 1 2β k ∥x k+1 -x∥ 2 M k+1 (1 + η k ) - 1 2 ∥ȳ k -y k ∥ 2 ≥ ϵ k + γσ k 2 ∥x k+1 -x∥ 2 .
(69)

Since δ ≤ 1, by dividing the above inequality with σ k , we have

1 2σ k ∥y k -ŷ∥ 2 - 1 2σ k ∥y k+1 -ŷ∥ 2 + 1 2τ k ∥x k -x∥ 2 M k - 1 2τ k ∥x k+1 -x∥ 2 M k+1 (1 + η k ) - 1 2σ k ∥ȳ k -y k ∥ 2 ≥ ϵ k σ k + γ 2 ∥x k+1 -x∥ 2 , (70)
where, we used τ k = β k σ k . By using (67), from the above we obtain that

1 2σ k ∥y k -ŷ∥ 2 - 1 2σ k ∥y k+1 -ŷ∥ 2 + 1 2τ k ∥x k -x∥ 2 M k - 1 2τ k ∥x k+1 -x∥ 2 M k+1 (1 + η k ) - 1 2σ k ∥ȳ k -y k ∥ 2 ≥ ϵ k σ k + γ 2C M ∥x k+1 -x∥ 2 M k+1 . (71)
It follows from the above inequality that

1 2σ k ∥y k -ŷ∥ 2 - 1 2σ k ∥y k+1 -ŷ∥ 2 + 1 2τ k ∥x k -x∥ 2 M k - 1 2σ k ∥ȳ k -y k ∥ 2 ≥ ϵ k σ k + 1 + (1 + η k )τ k γ/C M 2τ k (1 + η k ) ∥x k+1 -x∥ 2 M k+1 , 1 2σ k ∥y k -ŷ∥ 2 - 1 2σ k ∥y k+1 -ŷ∥ 2 + 1 2τ k ∥x k -x∥ 2 M k - 1 2σ k ∥ȳ k -y k ∥ 2 ≥ ϵ k σ k + τ k+1 (1 + τ k γ/C M ) τ k ∥x k+1 -x∥ 2 M k+1 2τ k+1 (1 + η k ) , (72)
For convenience, we set γ = γ/C M . From the update step of β k , it follows that

τ k+1 (1 + γτ k ) τ k ≥ τ k+1 min{C θ , (1 + γτ k )} τ k = σ k+1 σ k (73) Set B k := 1 2τ k ∥x k -x∥ 2 M k + 1 2σ k ∥y k -ŷ∥ 2 and Bk := B k Π k-1
i=1 (1+ηi) . From (72), we have:

σ k+1 σ k (1 + η k ) B k+1 + ϵ k σ k ≤ B k - 1 2σ k ∥ȳ k -y k ∥ 2 (74)
By dividing the above inequality by Π k-1 i=1 (1 + η i ) ≥ 1, we obtain

σ k+1 σ k Bk+1 + ϵ k σ k Π k-1 i=1 (1 + η i ) ≤ Bk - 1 2σ k Π k-1 i=1 (1 + η i ) ∥ȳ k -y k ∥ 2 (75)
By multiplying σ k on both sides, we have

σ k+1 Bk+1 + ϵ k Π k-1 i=1 (1 + η i ) ≤ σ k Bk - 1 2Π k-1 i=1 (1 + η i ) ∥ȳ k -y k ∥ 2 . ( 76
)
By Assumption 1, C = Π i∈N (1 + η i ) < +∞, we have

σ k+1 Bk+1 + ϵ k C ≤ σ k Bk - 1 2C ∥ȳ k y k ∥ 2 . ( 77
)
Summing up (77) from k = 1, • • • , N , we obtain

σ N +1 BN+1 + N k=1 ϵ k C ≤ σ 1 B1 - 1 2C N k=1 ∥ȳ k -y k ∥ 2 . ( 78
)
Since σ k is bounded by some σ for any k ∈ N, Bk is bounded from above. Since C = Π i∈N (1 + η i ) < +∞, B k is also bounded from above. So, y k is also bounded with lim k→∞ ∥ȳ k -y k ∥ 2 = 0. Thus, using the similar argument and notations in the proof B.2, we retrieve the same key inequality as the one in [START_REF] Malitsky | A first-order primal-dual algorithm with linesearch[END_REF]:

G( XN , Ȳ N ) ≤ C s N (σ 1 B 1 + θ 1 σ 1 P (x 0 )), ∥x N +1 -x∥ 2 M N +1 ≤ Cτ N +1 σ N +1 (σ 1 A 1 + θ 1 τ 1 P (x 0 )) = Cβ N +1 , (79)
Using the same argument from [START_REF] Malitsky | A first-order primal-dual algorithm with linesearch[END_REF], we know from B.1 that σ k is bounded by

µσ k = µ( -1+ √ (4δα)/β k +1 2 L
) where L = max{L, L K }. We claim that there exists a constant C β such that, β k = C β (1/k 2 ).

i If αδ/(β k ) ≤ 1, by σ k ≥ µσ k ≥ µσ, we have

β k+1 = β k min{C θ , 1 + γβ k σ k } ≤ β k min{C θ , 1 + µσδαγ} . (80)
In this case, β k decreases linearly. Thus, β k+1 ≤ C β /(k +1) 2 for k sufficiently large. ii If αδ/(β k ) ≥ 1, then σ k > µσ k > µ 2 L δα β k . Therefore, for k large enough, we have

β k+1 = β k min{C θ , 1 + γβ k σ k } ≤ β k min{C θ , 1 + µ √ δαγ 2 L √ β k } = β k 1 + µ √ δαγ 2 L √ β k . (

x∥ 2 M

 2 , and set prox τ g := prox I τ g ,

  where S k and Y k are matrices collecting the m most recent vectors in (12) as columns, D k := D(S ⊤ k Y k ) and L k := L(S ⊤ k Y k ) refer to the diagonal D(•) and the strict lower triangular L(•) part of the matrix S ⊤

  81)In this case, by inductionβ k ≤ C β k 2 for some constant C β > 0. From σ k > µσ k > µσ, we have s N = N k=1 σ k > N k=1 σ k > N k=1 O(k) ∼ N 2 since β k ≤ C β /k 2 fork sufficiently large. Then, we conclude the results.

⋆ We acknowledge funding by the ANR-DFG joint project TRINOM-DS under the number DFG OC150/5-1.

For example, the variable y k defined in[START_REF] Davis | Convergence rate analysis of primal-dual splitting schemes[END_REF] is not the dual variable in Algorithm 1.