Alternating Method of Successive Approximations

1. Solving Inverse Problems using Deep Learning

Deep learning (DL) techniques have witnessed remarkable progress in recent years, especially in the field of image processing, where they have consistently demonstrated significant improvements over traditional methods. Beyond image processing, DL and advanced optimization strategies have also found impactful applications in a range of complex domains such as wireless communication systems, quantum networking, and sensor network optimization, enabling solutions to previously intractable problems [1,2].

In the area of medical imaging, particularly for computed tomography (CT) and magnetic resonance imaging (MRI) reconstruction, one of the most successful DL-based strategies is known as algorithm unrolling or unfolding. This approach draws inspiration from classical iterative optimization algorithms, such as proximal gradient descent commonly used in variational methods. Instead of relying on handcrafted regularization terms designed through expert knowledge, unrolling methods integrate deep neural networks into the iterative framework. These networks are trained to learn optimal feature representations either from images or directly from sinogram data (the raw projection measurements), thereby enhancing the reconstruction process by capturing complex structures and noise patterns that traditional regularizers might miss.

More recently, dual-domain methods have gained attention, offering a further step forward by simultaneously utilizing complementary information from both the image domain and the sinogram (projection) domain. By bridging these two domains, these methods aim to improve reconstruction accuracy, resolution, and artifact suppression beyond what single-domain approaches can achieve.

Despite these exciting advances, DL-based reconstruction methods are not without limitations. One major challenge is their relatively weak theoretical foundation. Unlike classical optimization techniques that come with rigorous convergence and stability proofs, many DL models offer limited guarantees regarding their behavior during inference, particularly under distribution shifts or in clinical settings where robustness is critical. Moreover, these models are often memory-intensive, requiring substantial computational resources for training and inference, which can be a barrier to widespread adoption. Another persistent issue is the risk of overfitting, as the models may superficially imitate optimization procedures without truly capturing the underlying physics of the reconstruction problem.

Additionally, convolutional neural networks (CNNs), a cornerstone of deep learning, have demonstrated strong performance in several key tasks within medical imaging. They have been effectively employed for reconstruction from sparse-view and low-dose datasets, projection domain synthesis to enhance missing information, post-processing of preliminary reconstructions to refine image quality, and in the integration of learned priors into traditional iterative algorithms. In many cases, these DL-enhanced approaches have outperformed conventional analytical reconstruction methods, offering improvements in image quality, noise reduction, and diagnostic reliability.

In recent years, a new and increasingly influential class of deep learning (DL)-based methods, known as learnable optimization algorithms (LOAs), has been developed for image reconstruction tasks. Unlike earlier DL approaches, LOAs are designed with mathematical justifications and formal convergence guarantees, addressing some of the critical theoretical concerns traditionally associated with DL-based inverse problem solving. These methods have demonstrated significant advancements, particularly in areas such as MRI reconstruction [3].

LOAs are fundamentally rooted in the variational framework commonly used in classical inverse problem formulations [4]. In these methods, the regularization term—historically crafted manually based on prior knowledge—is instead parameterized using deep neural networks with learnable parameters. This transition allows the regularizer to capture much more complex and realistic image features, although it introduces challenges because the resulting objective function is often nonconvex and nonsmooth. Despite this complexity, LOAs aim to develop efficient, data-driven optimization schemes that ensure convergence towards a solution, typically by embedding algorithmic structures such as proximal gradient descent, half-quadratic splitting, or iterative shrinkage within the network design.

A distinctive characteristic of LOAs is that the network architecture closely mirrors the steps of an optimization algorithm: each layer corresponds to a logical iteration, and the entire deep network is highly structured according to the underlying optimization principles. Crucially, while the deep network’s parameters (such as proximal operators or penalty weights) are learned from data, the theoretical convergence properties—such as stability, robustness, and rates of convergence—are rigorously preserved, differentiating LOAs from more heuristic deep learning methods.

This innovative framework has been successfully employed to tackle practical challenges in medical imaging. For example, in sparse-view CT reconstruction, where the number of projection angles is limited to reduce radiation exposure, LOAs can adaptively learn regularization functions that suppress artifacts while preserving fine anatomical details. Similarly, in MRI reconstruction, LOAs have enabled faster imaging protocols by reconstructing high-quality images from significantly under-sampled k-space data.

Overall, learnable optimization algorithms represent a promising direction for combining the flexibility and expressive power of deep learning with the rigorous guarantees of traditional optimization theory, opening new opportunities for reliable and efficient image reconstruction across a wide range of applications.

Building upon prior work that established an optimal control framework for solving inverse problems with deep learnable regularizers [4], the current study extends this direction by integrating bilevel optimization with meta-learning strategies to further enhance adaptability and generalization across diverse tasks. By combining the dynamics of optimal control with meta-learned regularizers, our approach enables task-conditioned optimization that dynamically adjusts to different reconstruction scenarios. This synthesis offers a powerful strategy for improving both the theoretical rigor and practical performance of learning-based inverse problem solvers.

2. Optimal Control Viewpoint in Deep Learning

In this work, we propose an optimal control framework for solving inverse problems by minimizing a variational energy model with learned parameters. Specifically, we consider the following objective functional:

$$\mathcal{E}(\mathbf{u};\theta )=\mathcal{D}(\mathbf{u},\mathbf{b})+\lambda \mathcal{R}(\mathbf{u};\theta ),$$

(1)

where:

$\mathcal{D}(\mathbf{u},\mathbf{b})$ is the data fidelity term, enforcing consistency between the reconstructed signal $\mathbf{u}$ and the observed measurements $\mathbf{b}$;

$\mathcal{R}(\mathbf{u};\theta )$ is a regularization functional incorporating prior knowledge about $\mathbf{u}$, parameterized by a deep neural network (DNN) with learnable parameters $\theta$;

$\lambda >0$ balances the trade-off between data fidelity and regularization.

For linear inverse problems, the data fidelity is typically quadratic:

$$\mathcal{D}(\mathbf{u},\mathbf{b})={\displaystyle \frac{1}{2}}{|A\mathbf{u}-\mathbf{b}|}_2^2,$$

(2)

where A is a known measurement operator (e.g., a forward model for CT or MRI systems).

To minimize the variational energy $\mathcal{E}(\mathbf{u};\theta )$, we utilize the gradient flow dynamics, leading to the following ordinary differential equation (ODE) [5]:

$${\displaystyle \frac{d\mathbf{u}\left(t\right)}{dt}}=f(\mathbf{u}\left(t\right),\theta ):=-{\nabla }_{\mathbf{u}}\mathcal{E}(\mathbf{u}\left(t\right);\theta ),t\in [0,T],$$

(3)

where $T>0$ is a predefined terminal time. Here, $f\left(\mathbf{u}\right(t),\theta )$ represents the negative gradient of the energy functional at time t.

The goal is to learn the optimal parameters $\theta$ by casting the learning process as an optimal control problem (denoted as (P)), formulated as:

$$\underset{\theta \in \Theta }{min}J\left(\theta \right):=\Phi \left(\mathbf{u}\left(T\right)\right)+{\int }_0^{T}r(\mathbf{u}\left(t\right),\theta ),dt,$$

(4)

subject to the dynamical system constraint:

$$\left\{\begin{array}{cc}\dot{\mathbf{u}}\left(t\right)=f(\mathbf{u}\left(t\right),\theta ),\hfill & t\in [0,T],\mathbf{u}\left(0\right)={\mathbf{u}}^0,\hfill \end{array}\right.$$

(5)

where:

$\Phi \left(\mathbf{u}\right(T\left)\right)$ is the terminal cost, typically measuring the mismatch between the reconstructed final state and the ground truth;

$r\left(\mathbf{u}\right(t),\theta )$ is a running cost (often acting as a regularization or penalization term during the evolution);

$\Theta$ denotes the admissible set for control parameters;

${\mathbf{u}}^0$ represents the given initial condition.

The optimality conditions associated with this problem can be characterized by introducing a Lagrangian and deriving the corresponding adjoint state equations [6]. Define the Lagrangian $\mathcal{L}$:

$$\mathcal{L}(\mathbf{u},\mathbf{p},\theta )=\Phi \left(\mathbf{u}\left(T\right)\right)+{\int }_0^T\left(r(\mathbf{u}\left(t\right),\theta )+\langle \mathbf{p}\left(t\right),\dot{\mathbf{u}}\left(t\right)-f(\mathbf{u}\left(t\right),\theta )\rangle \right),dt,$$

(6)

where $\mathbf{p}\left(t\right)$ is the adjoint state.

Taking variations with respect to $\mathbf{u}$ and $\theta$ and setting them to zero leads to the following optimality system:

Minimizing the variational model (1) with the gradient flow method [7], we obtain the following ODE

$$\dot{\mathbf{u}}\left(t\right)=f(\mathbf{u}\left(t\right),\theta ):=-{\nabla }_{\mathbf{u}}\mathcal{E}(\mathbf{u}\left(t\right);\theta ),t\in [0,T].$$

(7)

We can cast the process of learning the undetermined parameters as an optimal control problem $\left(P\right)$ with control parameters $\theta$:

$$\begin{array}{cc}\hfill \underset{\theta \in \Theta }{min}& J\left(\theta \right):=\Phi \left(\mathbf{u}\left(T\right)\right)+{\int }_0^{T}r(\mathbf{u}\left(t\right),\theta )t,\left(P\right)\hfill \\ \hfill s.t.& \left\{\begin{array}{cc}\dot{\mathbf{u}}\left(t\right)=f(\mathbf{u}\left(t\right),\theta ),\hfill & 0\le t\le T,\hfill \\ \mathbf{u}\left(0\right)={\mathbf{u}}^0,\hfill \end{array}\right.\left(ODE\right)\hfill \end{array}$$

where $\Phi \left(\mathbf{u}\right(T\left)\right)$ is the terminal cost and $r\left(\mathbf{u}\right(t),\theta )$ is the running cost which usually plays role as the regularizer. Besides, $\Theta$ is the admissible set for the control parameters $\theta$ and ${\mathbf{u}}^0$ is the initial. These expressions enable the use of gradient-based optimization algorithms (e.g., stochastic gradient descent, Adam, or specialized adjoint-based methods) to update the network parameters $\theta$ efficiently while preserving the dynamics and structure of the underlying variational model.

Thus, our proposed framework unifies data-driven learning with the rigorous structure of variational optimization and control theory, providing a principled and theoretically grounded approach to solve challenging inverse problems such as sparse-view CT and accelerated MRI reconstruction.

3. Pontryagin’s Maximum Principle

First we give the definition of the Hamiltonian function $H:{\mathbb{R}}^n\times {\mathbb{R}}^n\times \Theta \to \mathbb{R}$ by

$$H(\mathbf{u};\mathbf{p};\theta )={\mathbf{p}}^{\top }f(\mathbf{u},\theta )-r(\mathbf{u},\theta ),$$

(8)

based on which we state the Pontryagin’s Maximum Principle (PMP) below

Theorem 1

(Pontryagin’s Maximum Principle, Informal Statement). If ${\theta }^{*}$ optimizes the optimal control problem(P), and ${\mathbf{u}}^{*}\left(t\right)$ is its corresponding state trajectory. Then, there exists an absolutely continuous co-state process ${\mathbf{p}}^{*}\left(t\right),t\in [0,T]$ such that

$$\left\{\begin{array}{c}{\dot{\mathbf{u}}}^{*}\left(t\right)={\nabla }_{\mathbf{p}}H({\mathbf{u}}^{*}\left(t\right);{\mathbf{p}}^{*}\left(t\right);{\theta }^{*}),{\mathbf{u}}^{*}\left(0\right)={\mathbf{u}}_0,\left(ODE\right)\hfill \\ {\dot{\mathbf{p}}}^{*}\left(t\right)=-{\nabla }_{\mathbf{u}}H({\mathbf{u}}^{*}\left(t\right);{\mathbf{p}}^{*}\left(t\right);{\theta }^{*}),{\mathbf{p}}^{*}\left(T\right)=-{\nabla }_{\mathbf{u}}\Phi {\left({\mathbf{u}}^{*}\left(T\right)\right)}^{\top },\left(ADJ\right)\hfill \\ H({\mathbf{u}}^{*}\left(t\right);{\mathbf{p}}^{*}\left(t\right);{\theta }^{*})=\underset{\theta }{max}H({\mathbf{u}}^{*}\left(t\right);{\mathbf{p}}^{*}\left(t\right);\theta ),\left(MAX\right)\hfill \end{array}\right.$$

are satisfied.

4. Algorithms of Successive Approximations

4.1. Basic Method of Successive Approximations

In light by the PMP, the basic Method of Successive Approximations (MSA) proposed in [8] are summarized in Algorithm 1 with an initial guess ${\theta }^0$ of the control parameter [4,9].

Algorithm 1 The Basic MSA [8] to Solve Inverse Problem

1: Initialize${\theta }^0$. 2: for$k=1,2,\cdots ,K$ (training iterations) do 3:    Solve ${\dot{\mathbf{u}}}^k\left(t\right)=f({\mathbf{u}}^k\left(t\right),{\theta }^{k-1})$, ${\mathbf{u}}^k\left(0\right)={\mathbf{u}}_0$; 4:    Solve ${\dot{\mathbf{p}}}^k\left(t\right)=-{\nabla }_{\mathbf{u}}H({\mathbf{u}}^k\left(t\right),{\mathbf{p}}^k\left(t\right),{\theta }^{k-1})$, ${\mathbf{p}}^k\left(T\right)=-\nabla \Phi {\left({\mathbf{u}}^k\left(T\right)\right)}^{\top }$; 5:    Set ${\theta }^k={error}_{\theta }{\int }_0^{T}H({\mathbf{u}}^k\left(t\right),{\mathbf{p}}^k\left(t\right),\theta )dt$; 6: end for 7: output:${\theta }^K$.

In practice, rather than exactly solving the maximization problem in Line 5, it is often more efficient to perform a gradient ascent step. Specifically, the control parameter $\theta$ is updated according to:

$${\theta }^k={\theta }^{k-1}+{\eta }_k{\nabla }_{\theta }{\int }_0^{T}H({\mathbf{u}}^k\left(t\right),{\mathbf{p}}^k\left(t\right),{\theta }^{k-1}),dt,$$

(9)

where ${\eta }_k>0$ denotes the learning rate (step size) at iteration k [4].

As proved in [8], with a gradient ascent (9) to maximize the Hamiltonian, the Algorithm 1 is equivalent to gradient descent with back-propagation (BP). The same as BP, one bottleneck of the basic MSA is the linear memory cost $\mathcal{O}\left(T\right)$ to cache all the intermediate states $\{({\mathbf{u}}^k\left(t\right),{\mathbf{p}}^k\left(t\right)):t\in [0,T]\}$. Despite its conceptual elegance, one major bottleneck of the basic MSA (similar to standard BP) is the linear memory complexity $\mathcal{O}\left(T\right)$ due to the need to store all intermediate states $({\mathbf{u}}^k\left(t\right),{\mathbf{p}}^k\left(t\right)):t\in [0,T]$ for accurate gradient computation. This memory burden can become significant for long time horizons or fine temporal discretizations, leading to scalability issues [4].

For practical implementation, the forward-backward system is typically discretized. Let us denote a discretization with N time steps, $t_n{n=0}^N$, where $t_n=n\Delta t$ and $\Delta t={\displaystyle \frac{T}{N}}$. Then, the forward propagation is:

$${\mathbf{u}}^{k}n+1={\mathbf{u}}_n^k+\Delta t,f({\mathbf{u}}_n^k,{\theta }^{k-1}),{\mathbf{u}}_0^k={\mathbf{u}}_0,$$

(10)

and the backward propagation becomes:

$${\mathbf{p}}_n^k={\mathbf{p}}_{n+1}^k+\Delta t,{\nabla }_{\mathbf{u}}H({\mathbf{u}}_n^k,{\mathbf{p}}_n^k,{\theta }^{k-1}),{\mathbf{p}}_N^k=-{\nabla }_{\mathbf{u}}\Phi {\left({\mathbf{u}}_N^k\right)}^{\top }.$$

(11)

The discrete Hamiltonian maximization step is then:

$${\theta }^k={\theta }^{k-1}+{\eta }_k\sum _{n=0}^{N-1}{\nabla }_{\theta }H({\mathbf{u}}_n^k,{\mathbf{p}}_n^k,{\theta }^{k-1})\Delta t.$$

(12)

These discrete updates highlight the need for efficient memory management strategies [4], such as checkpointing or reversible dynamics, to alleviate the $\mathcal{O}\left(T\right)$ memory overhead.

4.2. MSA with Augmented Reverse-State

As an alternative perspective, (9) can also be viewed as solving the following ODE backward in time

$${\theta }^k={\overline{\theta }}^k\left(0\right),where{\displaystyle \frac{d{\overline{\theta }}^k\left(t\right)}{dt}}={\eta }_k{\nabla }_{\theta }H({\mathbf{u}}^k\left(t\right),{\mathbf{p}}^k\left(t\right),{\overline{\theta }}^k\left(T\right)),with{\overline{\theta }}^k\left(T\right)={\theta }^{k-1}.$$

(13)

Suppose we integrate (13) backward w.r.t. time t, we get

$$\begin{array}{cc}\hfill {\overline{\theta }}^k\left(t\right)& ={\overline{\theta }}^k\left(T\right)+{\eta }_k{\int }_T^t{\nabla }_{\theta }H({\mathbf{u}}^k\left(t\right),{\mathbf{p}}^k\left(t\right),{\overline{\theta }}^k\left(T\right))dt\hfill \end{array}$$

(14a)

$$\begin{array}{cc}\hfill & ={\theta }^{k-1}+{\eta }_k{\int }_T^t{\nabla }_{\theta }H({\mathbf{u}}^k\left(t\right),{\mathbf{p}}^k\left(t\right),{\theta }^{k-1})dt(\mathrm{Because} \mathrm{we} \mathrm{set}{\overline{\theta }}^k\left(T\right)={\theta }^{k-1})\hfill \end{array}$$

(14b)

Then

$${\theta }^k={\overline{\theta }}^k\left(0\right)={\theta }^{k-1}+{\eta }_k{\int }_T^0{\nabla }_{\theta }H({\mathbf{u}}^k\left(t\right),{\mathbf{p}}^k\left(t\right),{\theta }^{k-1})dt$$

(15)

which is exactly the same as the gradient ascent of the Hamiltonian defined in (9).

From the definition of H, we know that ${\nabla }_{\theta }H(\mathbf{u},\mathbf{p},\theta )={\mathbf{p}}^{\top }{\nabla }_{\theta }f(\mathbf{u},\theta )$.

Instead of solving for co-state $\mathbf{p}$ and control $\theta$ separately, here we can solve for the augmented reverse-state $[\mathbf{p},\overline{\theta }]$ backward in time [4].

$${\theta }^k={\overline{\theta }}^k\left(0\right),where{\displaystyle \frac{d}{dt}}\left[\begin{array}{c}{\mathbf{p}}^k\left(t\right)\\ {\overline{\theta }}^k\left(t\right)\end{array}\right]=\left[\begin{array}{c}-{\nabla }_{\mathbf{u}}H({\mathbf{u}}^k\left(t\right),{\mathbf{p}}^k\left(t\right),{\overline{\theta }}^k\left(T\right))\\ {\eta }_k{\nabla }_{\theta }H({\mathbf{u}}^k\left(t\right),{\mathbf{p}}^k\left(t\right),{\overline{\theta }}^k\left(T\right))\end{array}\right],\left[\begin{array}{c}{\mathbf{p}}^k\left(T\right)\\ {\overline{\theta }}^k\left(T\right)\end{array}\right]=\left[\begin{array}{c}-\nabla \Phi {\left({\mathbf{u}}^k\left(T\right)\right)}^{\top }\\ {\theta }^{k-1}\end{array}\right],$$

(16)

which can be easily verified that its numerical result is identical to Algorithm 1 with Line 5 solved by (9). But the merit is that solving (16) only demands constant memory $\mathcal{O}\left(1\right)$ to cache the augmented reverse-state $[\mathbf{p},\overline{\theta }]$. However, in (16) we still need linear memory $\mathcal{O}\left(T\right)$ to store all intermediate forward state $\{{\mathbf{u}}^k\left(t\right):t\in [0,T]\}$. One way to tackle this issue is to augment it with the reverse state and solve for ${\mathbf{u}}^k\left(t\right)$ backward as well with backward initial ${\mathbf{u}}^k\left(T\right)$ computed from the forward pass [10], as shown in (17) [4].

$${\theta }^k={\overline{\theta }}^k\left(0\right),where{\displaystyle \frac{d}{dt}}\left[\begin{array}{c}{\mathbf{u}}^k\left(t\right)\\ {\mathbf{p}}^k\left(t\right)\\ {\overline{\theta }}^k\left(t\right)\end{array}\right]=\left[\begin{array}{c}{\nabla }_{\mathbf{p}}H({\mathbf{u}}^k\left(t\right),{\mathbf{p}}^k\left(t\right),{\overline{\theta }}^k\left(T\right))\\ -{\nabla }_{\mathbf{u}}H({\mathbf{u}}^k\left(t\right),{\mathbf{p}}^k\left(t\right),{\overline{\theta }}^k\left(T\right))\\ {\eta }_k{\nabla }_{\theta }H({\mathbf{u}}^k\left(t\right),{\mathbf{p}}^k\left(t\right),{\overline{\theta }}^k\left(T\right))\end{array}\right],\left[\begin{array}{c}{\mathbf{u}}^k\left(T\right)\\ {\mathbf{p}}^k\left(T\right)\\ {\overline{\theta }}^k\left(T\right)\end{array}\right]=\left[\begin{array}{c}{\mathbf{u}}^k\left(T\right)\\ -\nabla \Phi {\left({\mathbf{u}}^k\left(T\right)\right)}^{\top }\\ {\theta }^{k-1}\end{array}\right].$$

(17)

The whole process of (17) has constant memory cost $\mathcal{O}\left(1\right)$ which frees the space of $\{{\mathbf{u}}^k\left(t\right):t\in [0,T]\}$. But as a sacrifice we do trade the time for space as solving for (17) requires the re-computation of ${\mathbf{u}}^k\left(t\right)$.

If we replace the Line 4-5 in Algorithm 1 by (16) or (17), we can get an more memory-efficient algorithm with identical numerical result [4].

4.3. MSA with Backward Control Flow

As a careful observation in (13)-(17) on the right-hand side of the ODEs, all partial derivatives of H are evalauted on the initial state ${\overline{\theta }}^k\left(T\right)$. In this work, we further free this constraint and give more freedom to the dynamical system [4]. So instead of computing the partial derivative of H on the initial ${\overline{\theta }}^k\left(T\right)$, we compute it on the intermediate state ${\overline{\theta }}^k\left(t\right)$. As discussed in [11,12,13], the basic idea of the successive approximation methods is to find the optimal parameters from the guess by successive projections onto the manifold defined by the ODEs. So intuitively a better guess will contribute to better convergence performance and even better result. Our modification here gives a better guess to the optimal parameter ${\theta }^{*}$ for the Hamiltonian H at the intermediate time t, since ${\overline{\theta }}^k\left(t\right)$ is optimized backward in time so ${\overline{\theta }}^k\left(t\right)$ is usually a better estimation point than ${\overline{\theta }}^k\left(T\right)$. Along with (16), we summarize our proposed algorithm in Algorithm 2. The forward pass is to compute the trajectory of ${\mathbf{u}}^k\left(t\right)$, the backward pass is to compute the gradient flow for the co-state ${\mathbf{p}}^k\left(t\right)$ and control ${\overline{\theta }}^k\left(t\right)$. We would like to point out that here we did not extend (17) to the algorithm, because with changing to ${\overline{\theta }}^k\left(t\right)$, the trajectory of ${\mathbf{u}}^k\left(t\right)$ might flow to somewhere else, we cannot guarantee ${\mathbf{u}}^k\left(0\right)={\mathbf{u}}_0$ anymore. And this might cause the algorithm unstable. We will leave this problem to the future work [4].

Algorithm 2 The MSA with Reverse Augmented State and Control Flow to Solve Inverse Problem

1: Initialize${\theta }^0$. 2: for$k=1,2,\cdots ,K$ (training iterations) do 3:    Solve ${\dot{\mathbf{u}}}^k\left(t\right)=f({\mathbf{u}}^k\left(t\right),{\theta }^{k-1})$, ${\mathbf{u}}^k\left(0\right)={\mathbf{u}}_0$; 4:    Solve ${\displaystyle \frac{d}{dt}}\left[\begin{array}{c}{\mathbf{p}}^k\left(t\right){\overline{\theta }}^k\left(t\right)\end{array}\right]=\left[\begin{array}{c}-{\nabla }_{\mathbf{u}}H({\mathbf{u}}^k\left(t\right),{\mathbf{p}}^k\left(t\right),{\overline{\theta }}^k\left(t\right)){\eta }_k{\nabla }_{\theta }H({\mathbf{u}}^k\left(t\right),{\mathbf{p}}^k\left(t\right),{\overline{\theta }}^k\left(t\right))\end{array}\right],\left[\begin{array}{c}{\mathbf{p}}^k\left(T\right){\overline{\theta }}^k\left(T\right)\end{array}\right]=\left[\begin{array}{c}-\nabla \Phi {\left({\mathbf{u}}^k\left(T\right)\right)}^{\top }{\theta }^{k-1}\end{array}\right]$; 5:    Set ${\theta }^k={\overline{\theta }}^k\left(0\right)$. 6: end for 7: output:${\theta }^K$.

4.4. Time Discretization

We discretize $[0,T]$ to be $0,1,2,...,T$ and employ the explicit Euler method for the forward ODE and Verlet method [14] for the backward augmented ODE, then we have the discretized version of Algorithm 2, which is summarized in Algorithm 3. This algorithm was inspired by the learnable optimization algorithms for solving inverse problem [15,16], in the context of MRI reconstruction.

Algorithm 3 The Discretized Version of Algorithm 2 to Solve Inverse Problem [4]

1: Initialize${\theta }^0$. 2: for$k=1,2,\cdots ,K$ (training iterations) do 3:    Set ${\mathbf{u}}_0^k={\mathbf{u}}_0$; 4:    for $t=0,1,2,\cdots ,T-1$ do 5:      ${\mathbf{u}}_{t+1}^k={\mathbf{u}}_t^k+{\tau }_{k}f({\mathbf{u}}_t^k,{\theta }^{k-1})$; 6:    end for 7:    Set ${\mathbf{p}}_T^k=-\nabla \Phi {\left({\mathbf{u}}_T^k\right)}^{\top }$and${\overline{\theta }}_T^k={\theta }^{k-1}$; 8:    for $t=T-1,T-2,\cdots ,0$ do 9:      ${\overline{\theta }}_t^k={\overline{\theta }}_{t+1}^k+{\eta }_k{\nabla }_{\theta }H({\mathbf{u}}_t^k,{\mathbf{p}}_{t+1}^k,{\overline{\theta }}_{t+1}^k)$; 10:      ${\mathbf{p}}_t^k={\mathbf{p}}_{t+1}^k-{\nabla }_{\mathbf{u}}H({\mathbf{u}}_t^k,{\mathbf{p}}_{t+1}^k,{\overline{\theta }}_t^k)$; 11:    end for 12:    Set ${\theta }^{k+1}={\overline{\theta }}_0^{k+1}$; 13: end for 14: output:${\theta }^K$.

5. Design of the Regularizer in the Variational Model

In this study, we design the regularizer $\mathcal{R}$ to take the form $\mathcal{R}\left(\mathbf{u}\right)=\psi \left(\mathcal{G}\left(\mathbf{u}\right)\right)\in {\Re }^{+}$, where $\mathcal{G}:{\Re }^n\to {\Re }^n$ is a feature extraction operator achieved by deep neural network [4]. To enforce the sparsity, the function $\psi$ is defined as $\psi \left(\mathbf{u}\right)={\sum }_{i=1}^{n}logcosh\left({\mathbf{u}}_i\right)$, where $logcosh\left(x\right)$ is a twice differentiable function approximately equal to $\left|x\right|-log\left(2\right)$ for large x and to $x^2/2$ for small x. Throughout this work, we parameterize the feature extraction operator $\mathcal{G}$ as a vanilla l-layer convolutional neural network without bias separated by componentwise activation function as follows [4]:

$$\mathcal{G}\left(\mathbf{u}\right)={\mathbf{w}}_l*\sigma ...\sigma ({\mathbf{w}}_3*\sigma ({\mathbf{w}}_2*\sigma ({\mathbf{w}}_1*\mathbf{u}))),$$

(18)

where ${\left\{{\mathbf{w}}_k\right\}}_{k=1}^l$ denote the convolution weights and * denotes the convolution operation. Specifically, we parameterize the first convolution ${\mathbf{w}}_0$ to be d kernels of size $3\times 3$ and the last one ${\mathbf{w}}_l$ with 1 kernel of size $3\times 3\times d$. Besides, all hidden layers ${\left\{{\mathbf{w}}_k\right\}}_{k=2}^{l-1}$ correspond to convolutions with d kernels of size $3\times 3\times d$. Here, $\sigma$ represents a componentwise activation function which is twice differentiable. In this work we adopt the twice differentiable function sigmoid-weighted linear unit (SiLU) [17] as the activation.

6. Discussion

Recent advances in machine learning have significantly contributed to solving complex inverse problems and improving quantitative imaging techniques. In the context of MRI reconstruction, Bian et al. [18] introduced a self-supervised learning framework with model reinforcement, demonstrating enhanced accuracy and robustness for rapid $T_1$ mapping. Building upon these foundations, Bian et al. [19] proposed a diffusion modeling approach with domain-conditioned prior guidance, further accelerating both MRI and qMRI reconstructions while improving fidelity to the underlying physical models. From a theoretical perspective, Bian [4] also developed an optimal control-based bilevel optimization framework for inverse problems, providing a principled way to design deep learnable regularizers and control training dynamics.

Beyond medical imaging, in robotics and aerospace engineering, Gao et al. [20] proposed an autonomous multi-robot system for spacecraft servicing, showcasing the role of intelligent control and coordination in complex operational settings. Similarly, Yang et al. [21] explored the application of machine learning to assess digitalization capabilities in business finance, highlighting the versatility of data-driven methodologies across domains. Extending this direction, Gao et al. [22] addressed the adaptive detumbling of non-rigid satellites, integrating learning-based strategies to handle uncertainties in satellite dynamics. These advances collectively underscore the importance of combining domain knowledge, optimization theory, and machine learning techniques to address diverse challenges across imaging, finance, and autonomous systems.

7. Conclusions

This work is inspired by the work in [4] and recent advancements in multi-task MRI reconstruction using meta-learning approaches [23], we aim to design a scalable and principled learning framework that not only solves individual inverse problems efficiently but also adapts to varying data distributions and task-specific challenges.

In this work, we proposed a principled framework for solving inverse problems by casting deep learning-based reconstruction as an optimal control problem. Building on the variational model formulation, we introduced a bilevel optimization structure that integrates gradient flow dynamics with the Method of Successive Approximations (MSA) to train the network parameters in a theoretically grounded manner. Our analysis established connections between MSA and classical back-propagation, while highlighting the scalability challenges associated with memory growth. To address these limitations, we developed memory-efficient variants of MSA, including augmented reverse-state formulations and modified backward control flows, enabling constant memory cost while preserving convergence properties.

Through discrete-time implementations and the design of a learnable regularizer via deep convolutional networks, we demonstrated how this framework provides a flexible, scalable, and interpretable approach to tackle challenging inverse problems such as sparse-view CT and accelerated MRI reconstruction. Future directions include further improvements in algorithmic stability, integration with multi-physics models, and applications to broader classes of inverse problems beyond medical imaging.