[orcid=0009-0002-1984-1304]

[orcid=0000-0001-9296-0890]

Generalized Minimal Residual ,Tikhonov regularization ,Golub Kahan Bidiagonalization ,L-curve ,Generalized Cross-Validation

1 Introduction↩︎

X-ray computed tomography (CT) reconstructs an image of an object from X-ray measurements acquired at multiple projection angles. The reconstruction problem can be formulated as \[\label{eq:inv32prob} A\mathbf{x}\approx \mathbf{b}.\tag{1}\] Here \(\mathbf{x}\in \mathbb{R}^{n}\) is the unknown image, \(\mathbf{b}\in \mathbb{R}^{m}\) contains the measured projection data assumed to be contaminated by additive white Gaussian noise, and \(A \in \mathbb{R}^{m\times n}\) denotes the forward projection operator with no restrictions on its dimensions or rank.

The least squares problem for 1 \[\label{eq:LS} \min_{\mathbf{x}\in \mathbb{R}^n}\|A\mathbf{x}-\mathbf{b}\|_2,\tag{2}\] is solved by the normal equations, either \[\label{eq:over32NE} A^\top A\mathbf{x}=A^\top \mathbf{b},\tag{3}\] for \(m\geq n\), or \[\label{eq:under32NE} A A^\top \mathbf{y}= \mathbf{b}, \;\;\;\mathbf{x}=A^\top \mathbf{y},\tag{4}\] for the underdetermined case, \(m<n\). Both situations can be solved using standard iterative approaches, but for practical CT applications the dimensions of \(A\) are typically very large, making its explicit construction and storage computationally prohibitive. Consequently, iterative reconstruction algorithms are commonly implemented in a matrix-free framework, in which only matrix-vector products with \(A\) and its transpose are required. These operations are performed using numerical routines that approximate the underlying imaging physics via appropriate discretization schemes and the forward and back projection operators are often constructed using different discretizations and software implementations. As a result, the back projector \(B \in \mathbb{R}^{n\times m}\) is generally not the exact adjoint \(A^\top\) of the forward projector, that is, \(B \neq A^{\top}\). Then, the pair \(A,B\) is referred to as an unmatched projector pair [1], [2] and the normal equations are replaced with the generally unsymmetric unmatched normal equations, \[\label{eq:over32UN32NE} B A\mathbf{x}=B \mathbf{b},\tag{5}\] and \[\label{eq:under32UN32NE} A B \mathbf{y}= \mathbf{b}, \;\;\;\mathbf{x}=B \mathbf{y},\tag{6}\] respectively, for \(m \ge n\) and \(m<n\). If the range of \(A^\top\) equals the range of \(B\), then the solutions of the matched and unmatched equations are equivalent [3]. Here, we consider unmatched pairs that do not satisfy this range condition.

The use of unmatched projector pairs fundamentally alters the structure of the normal equations in 3 4 and classical iterative solvers, such as the Conjugate Gradient algorithm, do not necessarily yield reliable reconstructions [1], [4]. The AB- and BA-GMRES methods were introduced in [5] to solve general inverse problems with nonsymmetric operators, and several choices of the operator \(B\), including diagonal scaling strategies were designed to improve convergence. The AB- and BA-GMRES were adapted for CT reconstruction problems with unmatched projector pairs in [2]. There, GMRES is applied to the unmatched normal equations without requiring the explicit construction of the adjoint operator. Often exhibiting rapid initial convergence, these approaches remain susceptible to semiconvergence, in which early iterates approach a desirable approximation of the true solution before eventually becoming contaminated by propagated noise.

Incorporating Tikhonov regularization within the Krylov subspace updates reduces the issues of semiconvergence, providing the hybrid variants of AB- and BA-GMRES, [6]. Their effectiveness depends strongly on the properties of the projected problem generated by the Arnoldi process for which the Hessenberg matrix for the projected system governs both the regularization procedure and the quality of the reconstructed solution [7].

An alternative Krylov subspace framework is provided by the GKB process, originally developed by Golub and Kahan [8] for reducing a general matrix to bidiagonal form. Besides its importance in computing the singular value decomposition (SVD), GKB serves as the foundation of several widely used iterative methods, including LSQR [9] and LSMR [10]. Additional applications of GKB-based techniques can be found in [11]–[14]. Moreover, AB- and BA-GKB methods for reconstructing CT problems were presented in [15]. Unlike GMRES, which generates Hessenberg matrices, GKB produces lower bidiagonal matrices that are considerably sparser and retain important spectral information from the original problem. As a consequence, AB- and BA-GKB have been observed to exhibit milder semiconvergence and improved reconstruction quality compared with their GMRES counterparts [15].

Motivated by these observations, we develop hybrid AB- and BA-GKB methods for CT reconstruction with unmatched projector pairs. The bidiagonal structure of the projected space operator generated by GKB, assumed to be small, enables efficient construction of its SVD and facilitates efficient regularization parameter determination using standard methods. As examples, we consider L-curve and GCV for parameter estimation. Since regularization is applied only to a projected problem of relatively small dimension, the computational cost associated with the SVD and parameter-selection procedure is negligible compared with the matrix-vector operations required by the GKB iterations. Numerical results demonstrate that the proposed hybrid AB- and BA-GKB methods significantly reduce the effects of semiconvergence and can yield more accurate reconstructions than their GMRES-based counterparts.

1.1 Major Contributions↩︎

• Introduction of hybrid AB-GKB and BA-GKB algorithms for CT reconstruction with unmatched forward and backprojectors.

• Stabilization of the projected space operator using standard Tikhonov regularization, with GCV- and L-curve-based parameter selection using the SVD of the projected space operator.

• Numerical validation demonstrating that the proposed methods achieve reconstruction quality comparable to that of state-of-the-art hybrid GMRES approaches.

1.2 Organization↩︎

The organization of this paper is outlined as follows. 2 reviews the mathematical techniques used to solve inverse problems with unmatched pairs using GMRES, covering both classical and hybrid approaches. 3 revisits the GKB algorithm and introduces new \(AB\) and \(BA\) hybrid methodologies. 4 presents our numerical examples, while 5 concludes our study and highlights potential directions for future research.

2 AB- and BA- Iterative Methods for Unmatched Normal Equations↩︎

AB- and BA-GMRES methods for unmatched pairs were introduced in [2] and extended in the form as hybrid methods to include iterative regularization by [6]. Here we briefly review these algorithms for comparison with an alternative formulation and extension based on the GKB algorithm, [15].

AB-GMRES and BA-GMRES find at iteration \(k\) the least squares solution of \(\|AB\mathbf{y}-\mathbf{b}\|^2_2\), with update \(\mathbf{x}=B\mathbf{y}\), and of \(\|BA\mathbf{x}-\mathbf{b}\|^2_2\) for the unmatched normal equations 5 4 , respectively. Both algorithms use the same Krylov space \(\mathcal{K}_k(BA, B\mathbf{b})\) and construct an orthogonal basis for the approximate solution with least squares residual using an Arnoldi iteration. Let \(M\) be the coefficient matrix \(AB\) or \(BA\), then \(k\) applications of Arnoldi iteration yield \[\label{eq:Arn32iter} M \bar{Q}_k=\bar{Q}_{k+1}H_k,\tag{7}\] where \(\bar{Q}_{k+1}= [\bar{\mathbf{q}}_1,\bar{\mathbf{q}}_2,\dots,\bar{\mathbf{q}}_{k+1}]\;\) and the matrix \(H_{k+1,k}\) has upper Hessenberg structure, \[H_k= \begin{bmatrix}\label{upper32Hessenberg} h_{1,1}&h_{1,2} &\dots & h_{1,k}\\ h_{2,1}&h_{2,2} &\dots & h_{2,k}\\ &\ddots &\ddots &\vdots\\ &&h_{k,k-1}&h_{k,k}\\ &&&h_{k+1,k} \end{bmatrix}.\tag{8}\] Using this decomposition in 5 6 yields the projected problem \[\label{eq:proj} \mathbf{u}_k=\mathop{\mathrm{argmin}}_{\mathbf{u}} \Big\| H_k \mathbf{u}- \mathbf{e}_1 \|\mathbf{r}_0\|_2 \Big\|_2,\tag{9}\] where \(\mathbf{e}_1 =(1,0,\dots,0)^\top \in \mathbb{R}^{k+1}\).

For ill-conditioned problems, for which the underlying singular spectrum is broad, Krylov methods demonstrate semi-convergence. The initial iterations capture the dominant singular spectrum and the solution converges to the exact solution. But as the iteration proceeds, the Krylov subspace starts to include singular vectors, associated with smaller singular values, which are often contaminated by noise present in the data, and the solution departs from the exact solution. This problem is exacerbated for unmatched \(AB\) and \(BA\) operators for which the spectrum not only does not match the matched operator spectrum but is also unstable due to non-normality. Therefore, standard iterative methods applied to 9 typically require appropriate stopping rules to stop the iterations before the onset of semi-convergence. Then the number of iterations serves as a regularization parameter, [16]. Alternatively, one may seek a regularized solution on the projected space within the iteration, [17]. Applied to the AB- and BA-GMRES algorithms Tikhonov regularization on the projected space 9 yields the \(\lambda\)-dependent solution \[\label{eq:ABGMRES32Hyb} \mathbf{u}_{k,\lambda}=\mathop{\mathrm{argmin}}_{\mathbf{u}}\left \|\begin{bmatrix} H_k\\\lambda I \end{bmatrix} \mathbf{u}- \|\mathbf{r}_0\|_2\begin{bmatrix} \mathbf{e}_1\\ 0 \end{bmatrix} \right\|^2_2.\tag{10}\] An equivalent solution to 10 can be given using the solution of the normal equations \[\label{eq:normal32ABGMRES32Hyb} \mathbf{u}_{k,\lambda}= \|\mathbf{r}_0\|_2(H^\top_k H_k +\lambda^2I)^{-1} H^\top_k \mathbf{e}_1 .\tag{11}\] An option for the solution of 11 is to calculate the singular value decomposition of \(H_k\) and then use a standard parameter choice method to find \(\lambda\), [17].

The hybrid AB-GMRES and hybrid BA-GMRES Algorithms as detailed in [6] are summarized in [alg:ABGMRES,alg:BAGMRES].

The system matrix, \(H_k\), in [alg:ABGMRES,alg:BAGMRES] has potentially \((k^2+3k)/2\) non-zero elements and as \(k\) increases the computation of the SVD may become computationally expensive. In the next section, we review and extend an approach based on the use of the Golub Kahan Bidiagonalization (GKB) in place of GMRES, as introduced in [15].

3 Hybrid AB- and BA-GKB Algorithms↩︎

Following the work on the use of AB- and BA-GMRES for unmatched operators, AB- and BA-GKB were developed in [15]. Underlying the GKB algorithm is the decomposition \[\label{eq:GKB32dec} A Q_k=S_{k+1} C_k,\tag{12}\] where \(Q_{k}\) and \(S_{k+1}\) are column orthogonal matrices and \[\label{eq:lower32bid} C_{k} = \begin{bmatrix} \alpha_1&&&&\\ t_2&\alpha_2&&&\\ &\ddots&\ddots&&\\ &&&t_k&\alpha_k\\ &&&&t_{k+1} \end{bmatrix}\tag{13}\] is a lower bidiagonal matrix. For unmatched problems 5 6 , the GKB yields the projected problem, \(\min_{\mathbf{u}\in \mathbb{R}^k}\|C_k\mathbf{u}- t_1 \mathbf{e}_1\|_2^2\), corresponding to direct application of the GKB algorithm to the system matrix \(A\) with right, or left preconditioner \(B\), respectively. The solution is obtained as \(\mathbf{x}= Q_k\mathbf{u}_k\) in each case, but the AB- algorithm is initialized with \(\mathbf{b}/\|\mathbf{b}\|_2\) and the BA- algorithm with \(B\mathbf{b}/\|B\mathbf{b}\|_2\).

As for the algorithms based on the GMRES, an approach using the GKB is subject to semiconvergence of the iterations for ill-conditioned problems. Under the same framework with appropriate stopping conditions applied to GMRES and GKB algorithms for unmatched operators, [15] showed that the AB- and BA-GKB methods outperform their GMRES-based counterparts and are less susceptible to semiconvergence. To further improve reconstruction quality, we incorporate a regularization term into the AB- and BA-GKB methods yielding their hybrid counterparts based on the iterative solution of \[\mathbf{u}_{k,\lambda}=\mathop{\mathrm{argmin}}_{\mathbf{u}} \left\|\begin{bmatrix} \label{eq:ABGKB32Hyb} C_k \\ \lambda I \end{bmatrix} \mathbf{u}- t_1 \begin{bmatrix} \mathbf{e}_1\\ 0 \end{bmatrix} \right\|^2_2,\tag{14}\] equivalent to 10 but with bidiagonal \(C_k\) replacing upper Hessenberg \(H_k\), and where \(t_1=\|\mathbf{b}\|_2\) for the AB methods and \(\|B\mathbf{b}\|_2\) for the BA methods. The equivalent regularized solution to 14 based on the normal equations can be written as \[\label{eq:normal32ABGKB32Hyb} \mathbf{u}_{k,\lambda}= t_1(C^\top_k C_k +\lambda^2I)^{-1} C^\top_k \mathbf{e}_1 .\tag{15}\]

The hybrid AB-GKB and BA-GKB algorithms are summarized in [alg: Hybrid AB-GKB,alg: Hybrid BA-GKB].

3.1 Computational Costs of GMRES and GKB Algorithms↩︎

The computational costs of the non-hybrid GMRES and GKB algorithms were presented in [15] and are dominated by the costs of the matrix-vector operations. In particular, both the AB-GMRES and the AB-GKB algorithms require an extra matrix vector multiplication in the update of \(\mathbf{x}_k\), steps 14 and 15 in [alg:ABGMRES,alg: Hybrid AB-GKB], respectively. Thus, as compared to the respective BA algorithms the computational complexity increases by \(2mnk\). Supposing in each case, that the SVD of \(H_k\), respectively \(C_k\), is used to find the solution \(\mathbf{u}_{k,\lambda}\), then the difference between the GKB and GMRES algorithms is impacted by the structure of these matrices, the SVD for \(H_k\) being more expensive. As indicated in [15] the computational costs obtained are

\[\label{cost:ABGMRES} \text{Cost with AB-GMRES :} \mathcal{O}(6mnk+4mk^2+[(k^4+3k^2)/2]),\tag{16}\]

\[\label{cost:BAGMRES} \text{Cost with BA-GMRES :} \mathcal{O}(4mnk+4nk^2 + [(k^4+3k^2)/2] ),\tag{17}\]

\[\label{cost:ABGKB} \text{Cost with AB-GKB :} \mathcal{O}(10mnk+4mk^2+[2k^3]),\tag{18}\]

\[\label{cost:BAGKB} \text{Cost with BA-GKB :} \mathcal{O}(8mnk+4nk^2+ [2k^3]).\tag{19}\] While the terms in parentheses 16 17 and 18 19 , for the calculation of the SVD for \(H_k\) and \(C_k\), respectively, are small compared to large \(m\) and \(n\) for a limited number of iterations \(k\), it is evident that the SVD cost of \(C_k\) is approximately a factor of \(1/k\) smaller than that of \(H_k\).

Assuming that the algorithms stop with \(k<<[m,n]\), show that the GKB-based algorithms are respectively, five thirds and two times more computationally expensive that their GMRES counterparts for a fixed choice of \(k\). This was confirmed in [15] for the non-hybrid algorithms, where it was shown that the number of iterations to satisfaction of equivalent stopping rules was higher for the GKB algorithms and that their cost was higher. It is therefore of interest to investigate the performance of the comparative hybrid algorithms for the inversion of large scale problems. First, this requires a description of the approach used to find the regularization parameter \(\lambda\) at each iteration of each algorithm.

3.2 Selection of Regularization Parameters↩︎

The selection of \(\lambda\) in is important for the quality of the reconstruction. If the noise variance is known a priori, options to determine \(\lambda\) include the discrepancy principle or the unbiased predictive risk estimator. When the noise variance is unknown, methods for estimating \(\lambda\) include the L-curve and GCV. For consistency with real applications that often lack sufficient information about the noise variance, we consider the GCV, [18], and the L-curve, [19]. The motivation behind these choices is that our numerical experiments show that these methods effectively mitigate the impact of noise on reconstructions.

In our experiments we assume that \(k\) is small enough that we can use the SVD of \(C_k\), respectively \(H_k\), to evaluate the parameter-estimation techniques. We present the relevant formulae for the system with matrix \(C_k\) and right hand side \(t_1 \mathbf{e}_1\), noting that the formulae apply for \(H_k\) with right hand side \(\|\mathbf{r}_0\|_2\mathbf{e}_1\) in the same way. Now, let \(C_k = U \Sigma V\) denote the full SVD of \(C_k\), where \(U =[\mathbf{u}_1, \mathbf{u}_2, \dots,\mathbf{u}_k]\) and \(V=[\mathbf{v}_1, \mathbf{v}_2, \dots,\mathbf{v}_k]\) are orthogonal matrices and \(\Sigma\) is a (possibly rectangular) diagonal matrix containing the singular values \(\sigma_j\) on its main diagonal in descending order. The SVD filter factors take the form \(f_i(\lambda) = \sigma^2_i/(\sigma^2_i+\lambda^2)\), for \(i=1,2,\dots,k\), and the regularized solution \(\mathbf{u}_{k,\lambda}\), given in 15 , can be written as \[\label{eq:solnSVD}\mathbf{u}_{k,\lambda} = \sum^{k}_{i=1}t_1 f_i(\lambda) \frac{\mathbf{u}^\top_i\mathbf{e}_1}{\sigma_i}\mathbf{v}_i.\tag{20}\]

The GCV determines the minimum value of the GCV function, here given immediately in terms of the SVD as \[\label{eq:GCV32SVD} GCV(\lambda) = \frac{\sum^{k}_{i=1} t_1^2f^2_i(\lambda) (\mathbf{u}_i^\top\mathbf{e}_1)^2}{\sum^{k}_{i=1} f^2_i(\lambda)}.\tag{21}\] Even though the GCV is effective at finding a good estimate for \(\lambda\) for most problems, it has some difficulty when the GCV function is flat, [20]. This is because the noise in the data is small and the numerator in 21 is approximately constant for a large range of \(\lambda\). We consider, however, a reasonable amount of noise and the scenario of a flat GCV function is not a concern.

We also consider the L-curve in which \(\lambda\) is selected as the corner of the curve generated by plotting the log values of the following two quantities against each other, again in terms of the SVD \[\label{eq:SVD32Lcurve} \left( \sqrt{\sum^{k}_{i=1} t_1^2f_i(\lambda)^2 (\mathbf{u}^\top_i\mathbf{e}_1)^2} , \sqrt{\sum^{k}_{i=1} f_i(\lambda) t_1\frac{\mathbf{u}^\top_i\mathbf{e}_1}{\sigma_i}\mathbf{v}_i} \;\right).\tag{22}\] The corner in the L-curve indicates that the solution \(\mathbf{u}_k\) reaches a suitable balance between the noise remaining in \(\|\mathbf{u}_{k,\lambda}\|_2\) and the residual norm \(\|C_k\mathbf{u}_{k,\lambda}- t_1 \mathbf{e}_1\|_2\).

[eq:GCV SVD,eq:SVD Lcurve] facilitate efficient calculation of the regularization parameter over a range of values, based on the approximate singular values of the projected system matrix.

4 Numerical Experiments↩︎

In this section, we assess the performance of the proposed hybrid AB-GKB and BA-GKB algorithms for reconstructing CT problems with unmatched forward and backward projector pairs. We compare the reconstructions produced by these algorithms with those derived from hybrid AB-GMRES and BA-GMRES. For all simulations, we assume the availability of both the exact image \(\mathbf{x}^{\mathrm{true}}\) and the noise-free data \(\mathbf{b}^{\mathrm{true}}\). To evaluate the impact of noise, we introduce additive white Gaussian noise to \(\mathbf{b}^{\mathrm{true}}\). The Signal-to-Noise Ratio (SNR) is employed to quantify the noise present in \(\mathbf{b}\).

To evaluate the quality of the reconstruction, we record the relative reconstruction error (RRE) between the exact image \(\mathbf{x}^{\mathrm{true}}\) and the reconstructed image \(\mathbf{x}^{(k)}\) at iteration \(k\), as given by \[\label{eq:RE} \text{RRE}(\mathbf{x}_{k})=\frac{\|\mathbf{x}_{k}- \mathbf{x}^{\mathrm{true}}\|_2}{\|\mathbf{x}^{\mathrm{true}}\|_2}.\tag{23}\] In most practical problems, the exact image is not available; hence, using the RRE to stop the iterations at the point of semiconvergence is not possible. Here, we assume an accurate estimation of the noise level \(\|\mathbf{e}\|_2=\|\mathbf{b}^{\mathrm{true}}-\mathbf{b}\|_2\) is available and we terminate the iteration using the DP, [21], as \[\label{eq:DP} k_{DP}=\min_k{\|A\mathbf{x}_k-\mathbf{b}\|_2\leq \tau \|\mathbf{e}\|_2},\tag{24}\] where \(\tau\geq 1\) is a safety factor set to \(1\).

The Residual Norm Stagnation (RNS) occurs when the solver stops making meaningful progress. The iteration is terminated using RNS when \[\label{eq:RNS} k_{RNS}= \frac{ \bigr| \|A\mathbf{x}_{k-1}-\mathbf{b}\|_2 - \|A\mathbf{x}_k-\mathbf{b}\|_2 \bigr|}{\|A\mathbf{x}_{k-1}-\mathbf{b}\|_2} < \epsilon,\tag{25}\] where the tolerance \(\epsilon\) is user-defined to ensure convergence, and its appropriate choice may vary across algorithms. To facilitate a consistent comparison between the proposed methods and GMRES-based approaches, \(\epsilon\) is selected empirically for each simulation—based on numerical experiments (not reported)—to minimize the RRE for both hybrid AB-GMRES and BA-GMRES algorithms.

For each problem, we compare the results obtained using the DP in 24 with those from the RNS in 25 , with the RNS results always given in parentheses in the tables. In our implementations, we record the RRE and the corresponding iteration counts. Furthermore, we report average timings across \(5\) runs for \(150\) iterations to validate the theoretical computational discussion presented in 3.1. For each problem we also contrast the results using the GCV and L-curve for determination of the regularization parameter.

All forward and backward projectors used in this section are generated with the GPU-accelerated ASTRA Toolbox [22]. The operators \(A\) and \(B\) are implemented in a matrix-free manner, avoiding the explicit construction and storage of the corresponding matrices. Furthermore, the ASTRA Toolbox is extended to support the transpose operations \(A^\top\) and \(B^\top\), which are required for implementing the AB-GKB and BA-GKB algorithms.

All numerical results presented were computed using MATLAB Version 2026a [23] on a desktop equipped with an Intel(R) Core(TM) i5-14400F CPU, 16 GB of memory, and an NVIDIA GeForce RTX 5060 GPU. The software will be available on request to the author.

From 6 (a), it is immediately apparent that the hybrid BA-GMRES obtained using GCV is significantly affected by semiconvergence compared to all other methods. Furthermore, we can see that hybrid BA-GKB works well with GCV and yields excellent performance. The plots in 6 (b) illustrate that the L-curve offers more stable solutions that are comparable in RRE to those derived from GCV. However, the hybrid BA-GMRES with the L-curve still exhibits semiconvergence, but only for a higher number of iterations than when using the GCV. The results presented in column \(4\) of 1 support the effectiveness of DP against the RNS in terminating algorithms. Furthermore, we observe that while GCV and L-curve significantly affect the performance of the hybrid AB-GKB and AB-GMRES, this is not the case for their BA counterparts. We should emphasize that although \(\epsilon\) in the RNS is selected to minimize the GMRES methods, the GKB methods with DP can yield comparable RRE.

The execution times listed in column \(5\) of 1 show that the GMRES methods are faster than the GKB methods. Furthermore, the BA algorithms are faster than their AB counterparts. This validates the computational analysis presented in 3.1, .

The results reported in 2 demonstrate the capabilities of the hybrid BA-GKB algorithms to yield results comparable to those obtained with hybrid GMRES methods. We observe that the AB-GKB reaches the maximum number of iterations without convergence when using the GCV but that the BA-GKB has the same issue when using the L-curve. The recorded time in column \(5\) confirms the theoretical computational costs provided in 3.1, which show that the GMRES methods are cheaper than the GKB methods. For this larger underdetermined (\(m<<n\)) problem, it is clear that the BA algorithms are slower than the AB algorithms, generally requiring more iterations to convergence and the cost being dominated by the \(4nk^2>> 4mk^2\) term in in each case.

5 Conclusions↩︎

We introduced hybrid AB-GKB and BA-GKB algorithms for CT reconstruction with unmatched forward and backward projectors. The proposed methods were combined with GCV and L-curve parameter-selection strategies and evaluated using both the DP and RNS stopping criteria. Numerical results show that the proposed algorithms achieve reconstruction quality comparable to that of GMRES-based approaches. While GMRES remains computationally more efficient in terms of FLOPS, the proposed methods provide a viable alternative for solving unmatched-projector problems. No consistent superiority was observed between GCV and the L-curve methods, as their performance was problem-dependent. Moreover, the RNS stopping rule was generally satisfied faster than the DP stopping rule, but with a higher relative error. Future work will investigate additional regularization strategies and their integration with the proposed hybrid frameworks for both GMRES and GKB formulations.

References↩︎