<figure><img src="image_1.png" /></figure>

HAL
open science

# Interplay between preconditioning and regularization for linear ill-posed problems solved by conjugate gradient. Application to optical ow estimation

Ahmed Chabib, Jean-Francois Witz, Vincent Magnier, Pierre Gosselet

► To cite this version:

Ahmed Chabib, Jean-Francois Witz, Vincent Magnier, Pierre Gosselet. Interplay between preconditioning and regularization for linear ill-posed problems solved by conjugate gradient. Application to optical ow estimation. 2024. hal-04602735

HAL Id: hal-04602735

https://hal.science/hal-04602735

Preprint submitted on 6 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Conjugate gradient for ill-posed problems: regularization by preconditioning, preconditioning by regularization

Ahmed Chabib, Jean-François Witz, Vincent Magnier, Pierre Gosselet

► To cite this version:

Ahmed Chabib, Jean-Francois Witz, Vincent Magnier, Pierre Gosselet. Conjugate gradient for ill-posed problems: regularization by preconditioning, preconditioning by regularization. 2025. hal-04602735v2

HAL Id: hal-04602735

https://hal.science/hal-04602735v2

Preprint submitted on 10 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Interplay between preconditioning and regularization for linear ill-posed problems solved by conjugate gradient. Application to optical flow estimation

Ahmed Chabib, Jean-François Witz, Vincent Magnier, Pierre Gosselet

Univ. Lille, CNRS, Centrale Lille, UMR 9013 - LaMcube -
Laboratoire de Mécanique, Multiphysique, Multiéchelle,
F-59000 Lille, France

June 6, 2024

## Abstract.

This paper investigates the possibilities offered by combining regularization and preconditioning by the same symmetric positive semi-definite operator when solving ill-posed problems. We study the question of the stopping criterion, and the possibility offered by Ritz eigen elements for a posteriori filtering of the solution and tuning of Tikhonov's weight. The method is applied as the linear solver for an optical flow estimator and it is coupled with a subspace recycling strategy.

**Keywords.** regularization; preconditioning; conjugate gradient; Ritz values.

**License.** Published under CC-by license.

## 1 Introduction

Ill-posed systems of equations are ominous in mechanics. They are particularly present in identification problems, such as the boundary completion in elasticity [13, 7]. They also appear in methods involving some compact operator, like the Herglotz' transform to build solutions to the Helmholtz problems [14]. Beside issues of existence and uniqueness, ill-posed problems are characterized by the lack of stability between the cause and the effect, in other words small perturbations in the input potentially lead to large modifications of the output.

In this paper, we focus on discrete $n \times n$ linear symmetric positive semi-definite systems of the form $Ax = b$, allowing to analyze all properties in terms of the spectrum of $A$. Existence and uniqueness are linked to the null-space of $A$ (strictly zeros eigenvalues) whereas stability is associated with the accumulation of eigenvalues near zero. Indeed, a small contribution of $b$ in an eigendirection associated with a small eigenvalue of $A$ has a significant impact on $x$. Ill-posed problems thus result in poorly conditioned operators.

Solving such systems amounts to finding a satisfactory treatment to these small eigenvalues: truncation, shift, filtering. Truncation involves disregarding the problematic directions, for example using an eigenvalue decomposition $A = \sum_{i=1}^n \sigma_i u_i u_i^T$ (or more generally a singular value decomposition [9]) and only keeping the part of the matrix associated with eigenvalues larger than a given criterion $\epsilon$:

$$A^{-1} \approx \sum_{i=1}^{m} \frac{u_i u_i^T}{\sigma_i}, \quad \sigma_i > \epsilon > 0, \quad m < n. \qquad (1)$$
# Conjugate gradient for ill-posed problems: regularization by preconditioning, preconditioning by regularization

Ahmed Chabib, Jean-François Witz, Vincent Magnier, Pierre Gosselet
Univ. Lille, CNRS, Centrale Lille, UMR 9013 - LaMcube -
Laboratoire de Mécanique, Multiphysique, Multiéchelle,
F-59000 Lille, France

December 10, 2025

## Abstract

This paper investigates using the conjugate gradient iterative solver for ill-posed problems. We show that preconditioner and Tikhonov-regularization work in conjunction. In particular when they employ the same symmetric positive semi-definite operator, a powerful Ritz analysis allows one to estimate at negligible computational cost the solution for any Tikhonov's weight. This enhanced linear solver is applied to the boundary data completion problem and as the inner solver for the optical flow estimator. **Keywords:** regularization; preconditioning; conjugate gradient; Ritz values; L-curve; Picard plot.

## 1 Introduction

Ill-posed systems of equations are ubiquitous in mechanics. They are particularly present in identification problems, such as the boundary completion in elasticity [20, 13]. They also appear in methods involving some compact operator, like the Herglotz' transform to build solutions to the Helmholtz problems [21]. Beside issues of existence and uniqueness, ill-posed problems are characterized by the lack of stability between the cause and the effect, in other words small perturbations in the input potentially lead to large modifications of the output.

In this paper, we focus on discrete $n \times n$ linear(ized) symmetric positive semi-definite systems of the form $Ax = b$, allowing to analyze all properties in terms of the spectrum of $A$ which can be diagonalized as $A = \sum_{i=1}^{n} \sigma_i u_i u_i^T$ with $(u_i)$ an orthonormal basis and $(\sigma_i)$ non-negative eigenvalues in decreasing order.

Existence and uniqueness are linked to the null-space of $A$ (strictly zeros eigenvalues) whereas stability is associated with the accumulation of eigenvalues near zero. Indeed, a small contribution of $b$ in an eigendirection associated with a small eigenvalue of $A$ has a significant impact on $x = \sum \frac{u_i^T b}{\sigma_i} u_i$. Ill-posed problems thus result in poorly conditioned operators.

Solving such systems amounts to finding a satisfactory treatment to these small eigenvalues: truncation, shift, filtering.

* Truncation involves disregarding the problematic directions, and only keeping the part of the matrix associated with eigenvalues larger than a given criterion $\epsilon$:

$$A^{-1} \approx \sum_{i=1}^{m} \frac{u_i u_i^T}{\sigma_i}, \quad \sigma_i > \epsilon > 0 \text{ for } i \le m < n. \quad (1)$$

This method easily extends to general matrices using the singular value decomposition [15].

* Shift is generally achieved thanks to Tikhonov regularization [31], that can be written as:

$$A \approx A_{\lambda} = A + \lambda M. \quad (2)$$

The simplest choice is $M = I$ and in that case $\lambda > 0$ becomes the lower bound of the spectrum of $A_{\lambda}$. It is of course preferable to use a matrix $M$ with more physical sense, acting more locally on the small eigenvalues. Indeed, a “flat” regularization like $\lambda I$ may temper the contribution of the higher part of the spectrum. For instance one could use $\lambda(\sum_{i>m} u_i u_i^T)$, and the shift would behave like the truncation when $\lambda \to \infty$.

* Filtering tries either to “clean” the right-hand side from spurious excitation or to improve the solution after it was computed by enforcing some physical properties. For instance, smoothing can be used to recover
Note that this idea is approximately implemented by (non-preconditioned) iterative solvers since they tend to favor the upper part of the spectrum in the first iterations, so that one only needs to stop the solver early enough.

Shift is generally achieved thanks to Thikonov regularization [18], that can be written as, in its simplest form:

$$A \approx A_{\lambda} = A + \lambda I. \qquad (2)$$

In that case $\lambda > 0$ becomes the lower bound of the spectrum of $A_{\lambda}$. Often, a matrix with more physical sense, acting more locally on the small eigenvalues, is available instead of the identity.

Filtering tries to improve a solution after it was computed using another technique by enforcing some physical properties. For instance, smoothing can be used to recover regularity in a noisy solution.

All these techniques are often controlled by a parameter ($\epsilon$ for the truncation, $\lambda$ for the regularization, the stopping criterion of an iterative solver...) which needs to be tuned in order to find a balance between the information inside the original system and the information brought (or removed) by the treatment. When the accuracy of the data is known, Morozov's principle [15] provides an objective criterion for choosing the parameter: the correction introduced by the added information should not exceed the noise in the measurement.

When no such data is available, a compromise must be found. Picard's principle [10] compares the eigenvalues ($\sigma_i$) (sorted in decreasing order) and the decomposition of the right-hand side on the eigendirections ($u_i^T b$). While eigenvalues decrease less rapidly than their contributions to the right-hand side, the solution remains controlled. The L-curve [11] is a visual aid to find a balance. The solutions for various level of regularization are positioned in a frame ("norm of the residual", "norm of the solution"). In general large regularization leads to low norm of the solution but high error, whereas small regularization leads to lower level of error but large solutions (highly perturbed). Ideally, some corner exists which realizes a compromise between residual and oscillating solution.

In this paper, we attempt to combine these ideas within a sophisticated solver that offers many useful features: several stopping criteria, filtering of the solution and easy tuning of the regularization. It extends previous work in [7], by exploring the interplay between Tikhonov regularization and preconditioning. This work was initiated in the context of optical flow estimation, even though its scope is much broader. A key element in boosting the algorithm's capabilities is the ability to quickly solve the regularization operator, which is crucial in optical flow.

The paper is organized as follows. In Section 2, we recall the augmented preconditioned conjugate gradient algorithm and the computation of Ritz eigenelements, providing a first discussion on the effects of the preconditioner. In Section 3, we consider the case of regularized systems preconditioned by the regularization matrix. Section 4 provides assessments based on the optical flow estimation briefly explained in Appendix A. Section 5 concludes the paper.

## 2 Preconditioned Conjugate Gradient and Ritz elements

We use normal font for scalars, boldface lowercase for vectors and boldface uppercase for matrices. A collection of vectors ($x_j$) can be put in the matrix form $X_m = (x_0, ..., x_{m-1})$, the index $m$ thus corresponds to the number of columns of the matrix. We work in $\mathbb{R}^n$ even though the methods also apply for complex matrices and vectors.

Let $A$ be a symmetric definite positive matrix and $b$ be a vector. We search the solution to the system $Ax = b$. We use a conjugate gradient, preconditioned by the symmetric positive semi-definite matrix $M$, and augmented by the $n \times n_C$ full-rank matrix $C$ such that $ker(M) \subset Range(C)$.

At iteration $i$, we note $x_i$ the approximation and $r_i = b - Ax_i$ the residual. We introduce the augmented Krylov subspace $K_i(M^{-1}A, C, M^{-1}r_0)$ [5]:

$$K_i(M^{-1}A, C, M^{-1}r_0) = \operatorname{span} (M^{-1}r_0, \dots, (M^{-1}A)^{(i-1)}M^{-1}r_0) \oplus \operatorname{Range}(C) \quad (3)$$
regularity in a noisy solution. These approaches can be implemented by projectors: the former by right-
projection $\mathbf{x} = \mathbf{A}^{-1}\mathbf{P}_{r}\mathbf{b}$, and the latter by left projection $\mathbf{x} = \mathbf{P}_{l}\mathbf{A}^{-1}\mathbf{b}$. Note that the approaches are
equivalent when the preconditioner is orthogonal on an eigensubspace, $\mathbf{P} = \mathbf{P}^{T} = \sum_{i \le m} \mathbf{u}_{i}\mathbf{u}_{i}^{T}$ because
in that case $\mathbf{P}\mathbf{A} = \mathbf{A}\mathbf{P}$ and there is a direct correspondence between regular excitation and regular
solution. Another example of symmetric filtering is the use of coarse discretization in Galerkin methods.
For instance, in Digital Image Correlation [5], the coarser the discretization the better posed the problem
– but of course the less precise.

All these techniques are often controlled by a parameter (ε for the truncation, λ for the regularization,...)
which needs to be tuned in order to find a balance between the information inside the original system and
the information brought (or removed) by the treatment. When the accuracy of the data is known, Morozov's
principle [24] provides an objective criterion for choosing the parameter: the correction introduced by the added
information should not exceed the noise in the measurement.

When no such data is available, a compromise must be found. Picard's principle [16] compares the eigenvalues ($\sigma_i$) (sorted in decreasing order) and the decomposition of the right-hand side on the eigendirections ($\mathbf{u}_i^T \mathbf{b}$). While eigenvalues decrease less rapidly than their contributions to the right-hand side, the solution remains controlled. The L-curve [17] is a visual aid to find a balance. The solutions for various levels of regularization are positioned in a frame ("norm of the residual", "norm of the solution"). In general large regularization leads to low norm of the solution but high error, whereas small regularization leads to lower level of error but large solutions (highly perturbed). Ideally, some corner exists which realizes a trade-off between residual and oscillating solution.

Solving an ill-posed system with an iterative solver may seem counterintuitive, as poor preconditioning is
often considered a red flag for the use of such solvers. In fact, it can be viewed as an opportunity to implement
approximately the strategies mentioned above: because of their kinship with power iteration, iterative solvers
have the tendency to favor the higher part of the spectrum in the beginning, and the limited number of iterations
(often controlled by the convergence threshold) can be a way to stop iterations before searching too deeply in
the ill-conditioned portion of the spectrum. This paper focuses on the many opportunities opened by the
introduction of a preconditioner in the iteration.

Preconditioner $\mathbf{M}$ is often introduced in numerical methods courses as a cheap-to-compute approximation
of the inverse of the problem's operator ($\mathbf{M}^{-1} \approx \mathbf{A}^{-1}$) in the sense that the spectrum of $\mathbf{M}^{-1}\mathbf{A}$ should be as
concentrated as possible around a non-zero value (which can be scaled to 1). This can be roughly estimated by
the condition number of $\mathbf{M}^{-1}\mathbf{A}$, but more sophisticated studies are available [2]. However this does not apply
for ill-posed problems where the inverse of $\mathbf{A}$ can not be properly defined. Anyhow, the preconditioner may play
a major role as it can help tone the bad part of the spectrum down. It is well-known [29] that preconditioning
with SPD matrix $\mathbf{M}$ is equivalent to solving the system $\mathbf{L}^{-1}\mathbf{A}\mathbf{L}^{-T}\mathbf{y} = \mathbf{L}^{-1}\mathbf{b}$ where $\mathbf{L}\mathbf{L}^{T} = \mathbf{M}$ is the Cholesky
factorization, and $\mathbf{x} = \mathbf{L}^{-T}\mathbf{y}$. This system is governed by the eigenvalues of $\mathbf{L}^{-1}\mathbf{A}\mathbf{L}^{-T}$ or more practically by
the generalized eigenvalues ($\mu_i$) of ($\mathbf{A}, \mathbf{M}$): $\mathbf{A}\mathbf{v}_i = \mu_i\mathbf{M}\mathbf{v}_i$, with ($\mathbf{v}_i$) a $\mathbf{M}$-orthonormal basis. Using Rayleigh
quotient $R(\mathbf{v}) = (\mathbf{v}^T \mathbf{A} \mathbf{v}) / (\mathbf{v}^T \mathbf{M} \mathbf{v})$, one clearly sees that if the denominator is large enough in some directions
($\mathbf{v}$), the associated eigenvalues will be reduced and this part of the search space will be explored later in the
iterations. We see that the choice of the preconditioner obeys similar criteria as the choice of the regularization,
hence the use of the same notation $\mathbf{M}$.

In the case of well-posed problems, a precise solution is achievable whatever the preconditioner which mainly
impacts the number of iterations. In the case of ill-posed problem, preconditioning can be a tool to favor “better”
(more regular) directions. Since trying to achieve a too strict convergence may not be realistic, the sequence of
directions suggested by the preconditioner plays an important role in the definition of retained solution for a
relatively weak convergence threshold.

The use of preconditioned conjugate gradient seems not to be that frequent when solving ill-posed problems,
other Krylov solvers based of least-square (MinRes, CGNE, CGLS) [11, 32, 25, 7] or Landweber iteration [33]
are often preferred due to certain monotony properties. We believe that our study offer a set of important tools
to reconsider the potential of conjugate gradient for ill-posed problems.

## Main contributions

In this paper, we attempt to combine these ideas within a sophisticated preconditioned conjugate gradient
solver with several contributions. Note that we do not pretend to propose a new regularization method, but
we think that we bring new insight on the role of preconditioning and a particularly well-tuned solver where all
regularization tools can be deployed at marginal cost.

First we show that:

* The preconditioner can be leveraged to naturally regularize the solution.
Given an arbitrary initialization $\mathbf{x}_{00}$ and associated residual $\mathbf{r}_{00} = \mathbf{b} - \mathbf{A}\mathbf{x}_{00}$, the $i_{th}$ iteration can be defined as:

$$ \left\{ \begin{array}{ll} \text{find} & \mathbf{x}_i \in \mathbf{x}_{00} + \mathcal{K}_i(\mathbf{M}^{-1}\mathbf{A}, \mathbf{C}, \mathbf{M}^{-1}\mathbf{r}_{00}) \\ \text{such that} & \mathbf{r}_i \perp \mathcal{K}_i(\mathbf{M}^{-1}\mathbf{A}, \mathbf{C}, \mathbf{M}^{-1}\mathbf{r}_{00}) \end{array} \right. \qquad (4) $$

This iteration is achieved by Algorithm 1 where the augmentation is managed by the correction of the initialization (in order to obtain $\mathbf{x}_0$) and the projector $\mathbf{P}$ on $\ker(\mathbf{C}^T\mathbf{A})$, which together ensure that the residual remains orthogonal to $\text{Range}(\mathbf{C})$ [6].

<caption>**Algorithm 1 Augmented Conjugate Gradient**</caption>

$\mathbf{x}_{00}$ and $\mathbf{C}$ given
$\mathbf{P} = \mathbf{I} - \mathbf{C}(\mathbf{C}^T\mathbf{A}\mathbf{C})^{-1}\mathbf{C}^T\mathbf{A}$
$\mathbf{x}_0 = \mathbf{P}\mathbf{x}_{00} + \mathbf{C}(\mathbf{C}^T\mathbf{A}\mathbf{C})^{-1}\mathbf{C}^T\mathbf{b} = \mathbf{x}_{00} + \mathbf{C}(\mathbf{C}^T\mathbf{A}\mathbf{C})^{-1}\mathbf{C}^T\mathbf{r}_{00}$
$\mathbf{r}_0 = \mathbf{b} - \mathbf{A}\mathbf{x}_0 = \mathbf{P}^T\mathbf{r}_{00}$
$\mathbf{z}_0 = \mathbf{P}\mathbf{M}^{-1}\mathbf{r}_0, \mathbf{w}_0 = \mathbf{z}_0$
$\gamma_0 = (\mathbf{z}_0^T\mathbf{r}_0)$
**for** $i = 0, 1, \dots, m$ (convergence) **do**
&nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{q}_i = \mathbf{A}\mathbf{w}_i$
&nbsp;&nbsp;&nbsp;&nbsp;$\delta_i = (\mathbf{w}_i^T\mathbf{q}_i), \alpha_i = \delta_i^{-1}\gamma_i$
&nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{x}_{i+1} = \mathbf{x}_i + \mathbf{w}_i\alpha_i$
&nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{r}_{i+1} = \mathbf{r}_i - \mathbf{q}_i\alpha_i$
&nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{z}_{i+1} = \mathbf{P}\mathbf{M}^{-1}\mathbf{r}_{i+1}$
&nbsp;&nbsp;&nbsp;&nbsp;$\gamma_{i+1} = (\mathbf{z}_{i+1}^T\mathbf{r}_{i+1})$
&nbsp;&nbsp;&nbsp;&nbsp;$\beta_i = \gamma_i^{-1}\gamma_{i+1}$
&nbsp;&nbsp;&nbsp;&nbsp;$\mathbf{w}_{i+1} = \mathbf{z}_{i+1} + \mathbf{w}_i\beta_i$
**end for**

The algorithm builds two special basis of $\mathcal{K}_i(\mathbf{M}^{-1}\mathbf{A}, \mathbf{C}, \mathbf{M}^{-1}\mathbf{r}_{00})$, $\mathbf{Z}_i$ is $\mathbf{M}$-orthogonal whereas $\mathbf{W}_i$ is $\mathbf{A}$-orthogonal:

$$ \begin{aligned} \mathbf{Z}_i^T \mathbf{M} \mathbf{Z}_i &= \mathbf{Z}_i^T \mathbf{R}_i = \operatorname{diag}(\gamma_j)_{0 \le j < i} \\ \mathbf{W}_i^T \mathbf{A} \mathbf{W}_i &= \mathbf{W}_i^T \mathbf{Q}_i = \operatorname{diag}(\delta_j)_{0 \le j < i} \end{aligned} \qquad (5) $$

It is convenient to introduce the $\mathbf{M}$-normalized version of the $\mathbf{Z}_i$ basis:

$$ \hat{\mathbf{z}}_i = \frac{(-1)^i \mathbf{z}_i}{\sqrt{\gamma_i}} \quad \text{so that} \quad \hat{\mathbf{Z}}_i^T \mathbf{M} \hat{\mathbf{Z}}_i = \mathbf{I} \qquad (6) $$

$\hat{\mathbf{Z}}_i$ is in fact the basis that would have been obtained by the Arnoldi procedure [16], and we have:

$$ \begin{aligned} \hat{\mathbf{Z}}_i^T \mathbf{A} \hat{\mathbf{Z}}_i &= \mathbf{T}_i = \text{Tridiag}(\eta_{j-1}, \mu_j, \eta_j) \\ \text{with } \mu_0 = \frac{1}{\alpha_0}, \quad \mu_j = \frac{1}{\alpha_j} + \frac{\beta_{j-1}}{\alpha_{j-1}}, \quad \eta_j = \frac{\sqrt{\beta_j}}{\alpha_j} \end{aligned} \qquad (7) $$

We can diagonalize $\mathbf{T}_i = \Xi_i \Theta_i \Xi_i^T$ where $\Theta_i$ is the diagonal matrix of eigenvalues sorted in decreasing order and $\Xi_i$ the orthonormal matrix of eigenvectors.

The Ritz vectors are $\mathbf{V}_i = \hat{\mathbf{Z}}_i \Xi_i$, while $\Theta_i$ are the Ritz values of the system. They satisfy:

$$ \mathbf{V}_i^T \mathbf{M} \mathbf{V}_i = \mathbf{I} \quad \text{and} \quad \mathbf{V}_i^T \mathbf{A} \mathbf{V}_i = \mathbf{\Theta}_i. \qquad (8) $$

In order to mark the dependency of the Ritz vectors and values on the iteration $i$, they are denoted with an exponent (i): $\Theta_i = \operatorname{diag}(\theta_j^{(i)})_{1 \le j \le i}$ and $\mathbf{V}_i = (\mathbf{v}_1^{(i)}, \dots, \mathbf{v}_i^{(i)})$. As the number of iterations $i$ increases, the $(\theta_j^{(i)})_{1 \le j \le i}$ and $(\mathbf{v}_j^{(i)})_{1 \le j \le i}$ tend to approximate the generalized eigenvalues and eigenvectors of the couple $(\mathbf{A}, \mathbf{M})$ [12].
* The preconditioner provides appropriate norms to analyze the properties of the solution and plot ex-ploitable L-curves.

* In the spirit of [13] we show that the Ritz analysis may allow an interesting a posteriori filtering of the solution.

A key novelty lies in the cases where a Tikhonov regularization with an easily solvable structure is used, where we show that it is possible:

* To fully analyze the effect of the regularization on the original system in particular through Picard plots.

* To post-process at negligible cost an approximation of the solution for any other regularization weight $\lambda$.

When embedded in a nonlinear, process we show that:

* Acceleration is available.

* It is possible to postprocess approximations of the solutions for a predetermined family of weights ($\lambda_i$) at negligible cost.

For our assessments, we consider the Steklov-Poincaré approach to complete boundary data which has the advantage of being a linear problem, and the recovery of the optical flow between two images, which is nonlinear.

## Organization of the paper

The paper is organized as follows. In Section 2, we recall the augmented preconditioned conjugate gradient algorithm and the computation of Ritz eigenelements, we discuss and illustrate the effect of preconditioning and regularizing separately and the norms to analyze the evolution of the iteration. Particularly, in Section 3, we consider the case of regularized systems preconditioned by the regularization matrix where we can analyze in detail the effect of the regularization. These properties are assessed on the data completion problem in Section 4. Section 5 presents the adjustments required to handle nonlinear problems, in particular augmentation. Section 6 presents the assessments for estimating the optical flow. Section 7 concludes the paper.

## Notations

We use normal font for scalars, boldface lowercase for vectors and boldface uppercase for matrices. A collection of vectors ($\mathbf{x}_j$) can be put in the matrix form $\mathbf{X}_m = (\mathbf{x}_0, \dots, \mathbf{x}_{m-1})$, the index $m$ thus corresponds to the number of columns of the matrix. This work is presented in $\mathbb{R}^n$ even though the methods also apply for complex Hermitian matrices and vectors.

## 2 Preconditioned Conjugate Gradient and Ritz elements

Let $\mathbf{A}$ be a symmetric definite positive matrix and $\mathbf{b}$ be a vector. We search the solution to the system $\mathbf{A}\mathbf{x} = \mathbf{b}$. We use a conjugate gradient, preconditioned by the symmetric positive semi-definite matrix $\mathbf{M}$.

At iteration $i$, we note $\mathbf{x}_i$ the approximation and $\mathbf{r}_i = \mathbf{b} - \mathbf{A}\mathbf{x}_i$ the residual. We introduce the Krylov subspace $\mathcal{K}_i(\mathbf{M}^{-1}\mathbf{A}, \mathbf{M}^{-1}\mathbf{r}_0)$ [10]:

$$ \mathcal{K}_i(\mathbf{M}^{-1}\mathbf{A}, \mathbf{M}^{-1}\mathbf{r}_0) = \operatorname{span} (\mathbf{M}^{-1}\mathbf{r}_0, \dots, (\mathbf{M}^{-1}\mathbf{A})^{(i-1)}\mathbf{M}^{-1}\mathbf{r}_0) \quad (3) $$

Given an arbitrary initialization $\mathbf{x}_0$ and associated residual $\mathbf{r}_0 = \mathbf{b} - \mathbf{A}\mathbf{x}_0$, the $i_{th}$ iteration can be defined as:

$$ \left\{ \begin{array}{ll} \text{find} & \mathbf{x}_i \in \mathbf{x}_0 + \mathcal{K}_i(\mathbf{M}^{-1}\mathbf{A}, \mathbf{M}^{-1}\mathbf{r}_0) \\ \text{such that} & \mathbf{r}_i \perp \mathcal{K}_i(\mathbf{M}^{-1}\mathbf{A}, \mathbf{M}^{-1}\mathbf{r}_0) \end{array} \right. \quad (4) $$

This iteration is achieved by Algorithm 1.

The algorithm builds two special basis of $\mathcal{K}_i(\mathbf{M}^{-1}\mathbf{A}, \mathbf{M}^{-1}\mathbf{r}_0)$, $\mathbf{Z}_i$ is M-orthogonal whereas $\mathbf{W}_i$ is A-orthogonal:

$$ \mathbf{Z}_i^T \mathbf{M} \mathbf{Z}_i = \mathbf{Z}_i^T \mathbf{R}_i = \operatorname{diag}(\gamma_j)_{0 \le j < i} \quad (5) $$

$$ \mathbf{W}_i^T \mathbf{A} \mathbf{W}_i = \mathbf{W}_i^T \mathbf{Q}_i = \operatorname{diag}(\delta_j)_{0 \le j < i} $$

It is convenient to introduce the M-normalized version of the $\mathbf{Z}_i$ basis:

$$ \hat{\mathbf{z}}_i = \frac{(-1)^i \mathbf{z}_i}{\sqrt{\gamma_i}} \quad \text{so that} \quad \hat{\mathbf{Z}}_i^T \mathbf{M} \hat{\mathbf{Z}}_i = \mathbf{I} \quad (6) $$
## 2.1 Role of the preconditioner

It is often said that the preconditioner should be a good approximation of the inverse of A, in the sense that the spectrum of M⁻¹A should be as concentrated as possible around a non-zero value (which can be scaled to 1). This can be roughly estimated by the condition number of M⁻¹A, but more sophisticated studies are available [2]. It is important to note the proximity between the conjugate gradient algorithm to solve linear systems and the Lanczos procedure to compute eigenvalues [17], as was made explicit by the Ritz analysis in this section’s introduction. It is also useful to see that the higher part of the spectrum is explored in priority due to the repeated power in the construction of Krylov subspace.

In the case of poorly-conditioned systems, the preconditioner can play a regularization role, as was explored in [7]. Schematically, for a direction d, what matters is the ratio (d^TAd/d^TMd). If the preconditioner measures the irregularity of a field, it penalizes the highly oscillating directions and delays their exploration.

Also, the preconditioner can be viewed as providing a physic-based alternative to the simple Euclidean orthogonality. It thus defines “natural” norms which are useful in the analysis of the iterations, as discussed in next subsection.

## 2.2 Stopping criteria

Conjugate gradient gives valuable pieces of information at no cost, but in specific norms, in the course of the iterations. First, we have error estimators [1]:

$$ 
\begin{aligned}
\|\mathbf{r}_i\|_{M^{-1}}^2 &= \gamma_i \\
\|\mathbf{x}_{i+1} - \mathbf{x}\|_{A}^2 &= \|\mathbf{x}_i - \mathbf{x}\|_{A}^2 - \gamma_i^2 \delta_i^{-1}
\end{aligned}
\qquad (9) $$

of course the difficulty for the second identity is that $\|\mathbf{x}_0 - \mathbf{x}\|_{A}^2$ is unknown. We also have measurement of the norm of the correction brought by iterations [7]:

$$ 
\begin{aligned}
\|\mathbf{x}_{i+1} - \mathbf{x}_0\|_{M}^2 &= \|\mathbf{x}_i - \mathbf{x}_0\|_{M}^2 + \alpha_i^2 \|\mathbf{w}_i\|_{M}^2 + 2\alpha_i (\mathbf{w}_i^T \mathbf{M} (\mathbf{x}_i - \mathbf{x}_0)) \\
&\quad \text{with } \begin{cases} \|\mathbf{w}_{i+1}\|_{M}^2 = \gamma_i + \beta_i^2 \|\mathbf{w}_i\|_{M}^2, & \|\mathbf{w}_0\|_{M}^2 = \gamma_0, \\ (\mathbf{w}_{i+1}^T \mathbf{M} (\mathbf{x}_{i+1} - \mathbf{x}_0)) = -\beta_i ((\mathbf{w}_i^T \mathbf{M} (\mathbf{x}_i - \mathbf{x}_0)) + \alpha_i \|\mathbf{w}_i\|_{M}^2). \end{cases}
\end{aligned}
\qquad (10) $$

Finally, we have an estimator on the preconditioned operator:

$$ 
\begin{aligned}
\|\mathbf{T}_0\|_F^2 &= \mu_0^2, \\
\|\mathbf{T}_{i+1}\|_F^2 &= \|\mathbf{T}_i\|_F^2 + \mu_i^2 + \eta_i^2 + \eta_{i-1}^2 \rightarrow \|\mathbf{M}^{-1}\mathbf{A}\|_F^2,
\end{aligned}
\qquad (11) $$

where index F stands for the Frobenius norm, $\|\mathbf{M}^{-1}\mathbf{A}\|_F^2$ is the sum of the squares of the generalized eigenvalues of (A, M).

We can then devise costless stopping criteria:

$$ 
\begin{aligned}
\|\mathbf{r}_i\|_{M^{-1}}^2 &< \varepsilon \|\mathbf{r}_0\|_{M^{-1}}, \\
\|\mathbf{r}_i\|_{M^{-1}}^2 &< \varepsilon \|\mathbf{T}_i\|_F \|\mathbf{x}_i - \mathbf{x}_0\|_{M}.
\end{aligned}
\qquad (12) $$

The first one is very classical, but it is risky in the sense that it may be too strict is the initialization or the augmentation were well-chosen ($\|\mathbf{r}_0\|_{M^{-1}}$ is already small). The second one is inspired from the Scipy implementation of MinRes with a more adapted choice of norms, we are sorry not to know whom to attribute it to. It is interesting in the sense that it balances the reduction of error and the increase of the norm of the solution, which is always a dilemma when solving ill-posed problems. It is often interesting to combine the criteria, add stagnation detection, and to also use safeguards in absolute value in case of too good initialization and augmentation.
<caption>**Algorithm 1** Preconditioned Conjugate Gradient</caption>

$\mathbf{r}_0 = \mathbf{b} - \mathbf{A}\mathbf{x}_0 = \mathbf{P}^T\mathbf{r}_0$
$\mathbf{z}_0 = \mathbf{P}\mathbf{M}^{-1}\mathbf{r}_0, \mathbf{w}_0 = \mathbf{z}_0$
$\gamma_0 = (\mathbf{z}_0^T\mathbf{r}_0)$
**for** $i = 0, 1, \dots, m$ (convergence) **do**
&nbsp;&nbsp;&nbsp;&nbsp;**q**$_i = \mathbf{A}\mathbf{w}_i$
&nbsp;&nbsp;&nbsp;&nbsp;$\delta_i = (\mathbf{w}_i^T\mathbf{q}_i), \alpha_i = \delta_i^{-1}\gamma_i$
&nbsp;&nbsp;&nbsp;&nbsp;**x**$_{i+1} = \mathbf{x}_i + \mathbf{w}_i\alpha_i$
&nbsp;&nbsp;&nbsp;&nbsp;**r**$_{i+1} = \mathbf{r}_i - \mathbf{q}_i\alpha_i$
&nbsp;&nbsp;&nbsp;&nbsp;**z**$_{i+1} = \mathbf{P}\mathbf{M}^{-1}\mathbf{r}_i$
&nbsp;&nbsp;&nbsp;&nbsp;$\gamma_{i+1} = (\mathbf{z}_{i+1}^T\mathbf{r}_i)$
&nbsp;&nbsp;&nbsp;&nbsp;$\beta_i = \gamma_i^{-1}\gamma_{i+1}$
&nbsp;&nbsp;&nbsp;&nbsp;**w**$_{i+1} = \mathbf{z}_{i+1} + \mathbf{w}_i\beta_i$
**end for**

$\hat{\mathbf{Z}}_i$ is in fact the basis that would have been obtained by the Arnoldi procedure [28], and we have:

$$
\begin{gathered}
\hat{\mathbf{Z}}_i^T \mathbf{A} \hat{\mathbf{Z}}_i = \mathbf{T}_i = \text{Tridiag}(\eta_{j-1}, \mu_j, \eta_j) \\
\text{with } \mu_0 = \frac{1}{\alpha_0}, \quad \mu_j = \frac{1}{\alpha_j} + \frac{\beta_{j-1}}{\alpha_{j-1}}, \quad \eta_j = \frac{\sqrt{\beta_j}}{\alpha_j}
\end{gathered}
\tag{7}
$$

We can diagonalize $\mathbf{T}_i = \Xi_i \mathbf{\Theta}_i \Xi_i^T$ where $\mathbf{\Theta}_i$ is the diagonal matrix of eigenvalues sorted in decreasing order and $\Xi_i$ the orthonormal matrix of eigenvectors.

The Ritz vectors are $\mathbf{V}_i = \hat{\mathbf{Z}}_i \Xi_i$, while $\mathbf{\Theta}_i$ are the Ritz values of the system. They satisfy:

$$
\mathbf{V}_i^T \mathbf{M} \mathbf{V}_i = \mathbf{I} \quad \text{and} \quad \mathbf{V}_i^T \mathbf{A} \mathbf{V}_i = \mathbf{\Theta}_i. \tag{8}
$$

In order to mark the dependency of the Ritz vectors and values on the iteration $i$, they are denoted with an exponent (i): $\mathbf{\Theta}_i = \text{diag}(\theta_j^{(i)})_{1\le j\le i}$ and $\mathbf{V}_i = (\mathbf{v}_1^{(i)}, \dots, \mathbf{v}_i^{(i)})$. As the number of iterations $i$ increases, the $(\theta_j^{(i)})_{1\le j\le i}$ and $(\mathbf{v}_j^{(i)})_{1\le j\le i}$ tend to approximate the generalized eigenvalues and eigenvectors of the couple $(\mathbf{A}, \mathbf{M})$ [19].

## 2.1 Appropriate solution and error norms and stopping criteria

Conjugate gradient gives valuable pieces of information at no cost in the course of the iterations, but in specific
norms. Indeed, the preconditioner can be viewed as providing a physic-based alternative to the simple Euclidean
norm and inner product.

First, we have error estimators [18, 1]:

$$
\begin{align}
\|\mathbf{r}_i\|_{\mathbf{M}^{-1}}^2 &= \gamma_i \nonumber \\
\|\mathbf{x}_{i+1} - \mathbf{x}\|_{\mathbf{A}}^2 &= \|\mathbf{x}_i - \mathbf{x}\|_{\mathbf{A}}^2 - \gamma_i^2 \delta_i \tag{9}
\end{align}
$$

of course the difficulty for the second identity is that $\|\mathbf{x}_0 - \mathbf{x}\|_{\mathbf{A}}^2$ is unknown. Note that in [23], a strategy is proposed to overcome this problem and to fully estimate $\|\mathbf{x}_{i+1} - \mathbf{x}\|_{\mathbf{A}}^2$.

We also have measurement of the norm of the correction brought by iterations [13]:

$$
\begin{align}
\|\mathbf{x}_{i+1} - \mathbf{x}_0\|_{\mathbf{M}}^2 &= \|\mathbf{x}_i - \mathbf{x}_0\|_{\mathbf{M}}^2 + \alpha_i^2 \|\mathbf{w}_i\|_{\mathbf{M}}^2 + 2\alpha_i (\mathbf{w}_i^T \mathbf{M} (\mathbf{x}_i - \mathbf{x}_0)) \nonumber \\
&\quad \text{with } \begin{cases} \|\mathbf{w}_{i+1}\|_{\mathbf{M}}^2 = \gamma_i + \beta_i^2 \|\mathbf{w}_i\|_{\mathbf{M}}^2, & \|\mathbf{w}_0\|_{\mathbf{M}}^2 = \gamma_0, \\ (\mathbf{w}_{i+1}^T \mathbf{M} (\mathbf{x}_{i+1} - \mathbf{x}_0)) = -\beta_i ((\mathbf{w}_i^T \mathbf{M} (\mathbf{x}_i - \mathbf{x}_0)) + \alpha_i \|\mathbf{w}_i\|_{\mathbf{M}}^2). \end{cases} \tag{10}
\end{align}
$$

Finally, we have an estimator on the preconditioned operator:

$$
\begin{align}
\| \mathbf{T}_0 \|_F^2 &= \mu_0^2, \nonumber \\
\| \mathbf{T}_{i+1} \|_F^2 &= \| \mathbf{T}_i \|_F^2 + \mu_i^2 + \eta_i^2 + \eta_{i-1}^2 \rightarrow \| \mathbf{M}^{-\frac{1}{2}} \mathbf{A} \mathbf{M}^{-\frac{1}{2}} \|_F^2, \tag{11}
\end{align}
$$

where index *F* stands for the Frobenius norm, $\|\mathbf{M}^{-\frac{1}{2}}\mathbf{A}\mathbf{M}^{-\frac{1}{2}}\|_F^2$ is the sum of the squares of the generalized eigenvalues of $(\mathbf{A}, \mathbf{M})$.

We can then devise costless (without extra computation) stopping criteria:

$$
\|\mathbf{r}_i\|_{\mathbf{M}^{-1}} < \varepsilon \|\mathbf{r}_0\|_{\mathbf{M}^{-1}}, \qquad (12a)
$$
## 2.3 A posteriori filtrering

In the case of a poorly conditioned system, the reduction of the error can be obtained at the price of an explosion of the norm of the solution. This is well explained by Picard analysis: the phenomenon occurs when the eigenvalues of the operator decrease faster than the contribution of the right-hand side in the associated direction. It can also be visualized on a L-curve, in the positive quarter of a frame of the form ($\|\mathbf{r}_i\|, \|\mathbf{x}_{i+1}\|$): the curve starts in the bottom right corner (large error, small norm) with a fast decay of the error, and finishes in the top left corner (reduced error, large norm).

As shown earlier, conjugate gradient provides natural norms to evaluate the error and the norm or the solution: $\|\mathbf{x}_i - \mathbf{x}\|_A$ and $\|\mathbf{x}_i - \mathbf{x}_0\|_M$. With this choice of norms, the curve is always oriented toward the upper-left corner: at each iteration, the norm of the error decreases and the norm of the solution increases.

Ritz elements offer a convenient way to filter the solution. Assuming $m$ iterations were conducted, we can process the basis $\mathbf{V}_m$ and the values $\Theta_m$. We can decompose the right-hand side on the Ritz basis $r_j^{(m)} = \mathbf{v}_j^{(m)T} \mathbf{r}_0$, and define:

$$ \text{for } i \le m, \quad \tilde{\mathbf{x}}_i^{(m)} = \mathbf{x}_0 + \sum_{j=1}^{i} \frac{r_j^{(m)}}{\theta_j^{(m)}} \mathbf{v}_j^{(m)}. \qquad (13) $$

We have:

$$ \|\tilde{\mathbf{x}}_i^{(m)} - \mathbf{x}\|_A^2 = \|\tilde{\mathbf{x}}_i^{(m)} - \mathbf{x}_0\|_A^2 - \sum_{j=1}^{i} \frac{(r_j^{(m)})^2}{\theta_j^{(m)}}, \qquad (14) $$

$$ \|\tilde{\mathbf{x}}_i^{(m)} - \mathbf{x}_0\|_M^2 = \sum_{j=1}^{i} \frac{(r_j^{(m)})^2}{(\theta_j^{(m)})^2}, $$

and of course:

$$ \|\tilde{\mathbf{x}}_i^{(m)} - \tilde{\mathbf{x}}_{i-1}^{(m)}\|_A^2 = \frac{(r_i^{(m)})^2}{\theta_i^{(m)}} \quad \text{and} \quad \|\tilde{\mathbf{x}}_i^{(m)} - \tilde{\mathbf{x}}_{i-1}^{(m)}\|_M^2 = \frac{(r_i^{(m)})^2}{(\theta_i^{(m)})^2}. \qquad (15) $$

Since the $(\theta_j^{(m)})$ are sorted in decreasing order, we see that the error of $(\tilde{\mathbf{x}}_i^{(m)})$ tends to decrease slower than its norm tends to increase. The L-curve for $(\tilde{\mathbf{x}}_i^{(m)})_i$ is then convex and the corner may be easier to define. The slope of the L-curve between the point $i-1$ and $i$ is $-(\theta_i^{(m)})^{-1}$. A possibility is to define the corner as the point which maximizes the variation of slope: $i = \arg\max_j((\theta_{j+1}^{(m)})^{-1} - (\theta_j^{(m)})^{-1})$.

Ritz' elements also make it possible to use Picard's theory and stop the construction of $\tilde{\mathbf{x}}_i^{(m)}$ when the contribution $r_j^{(m)}$ starts to decrease less fast than $\theta_j^{(m)}$. This criterion has the advantage to take into account the properties of the right-hand side.

# 3 Preconditioning by regularization

In general, preconditioning is important because it can speed up the convergence of iterative solvers. As evoked in Section 2, it becomes crucial for ill-posed problem for the physical information it provides and also because the resolution will be stopped “early”. Since a too strict convergence criterion can not be attained in a reasonable amount of time, two distinct preconditioners lead to different history of resolution and to qualitatively very different solutions.

We are interested in Tikhonov-regularized systems of the form:

$$ \underbrace{(\mathbf{A} + \lambda\mathbf{M})\mathbf{x}_{\lambda}}_{\mathbf{A}_{\lambda}} = \underbrace{\mathbf{b}_{\mathbf{A}} + \lambda\mathbf{b}_{\mathbf{M}}}_{\mathbf{b}_{\lambda}} \qquad (16) $$
$$ \|\mathbf{r}_i\|_{\mathbf{M}^{-1}} < \varepsilon \|\mathbf{T}_i\|_F \|\mathbf{x}_i - \mathbf{x}_0\|_{\mathbf{M}}, \qquad (12b) $$

$$ \gamma_i^2 \delta_i^{-1} < \varepsilon^2. \qquad (12c) $$

The first one is very classical, but it is risky in the sense that it may be too strict if the initialization was well-chosen ($\|\mathbf{r}_0\|_{\mathbf{M}^{-1}}$ is already small). The second one is inspired from the Scipy implementation of MinRes with a more adapted choice of norms, we were unable to trace the original source of this idea. This criterion is useful because it balances the error reduction against the growth of the solution norm, a central dilemma when solving ill-posed problems. The third criterion corresponds to the iteration bringing negligible correction, what is often called stagnation, in general it must be checked for a sequence of iterations before actually stopping. It is often interesting to combine the criteria, add stagnation detection, and to also use safeguards in absolute value in case of too good initialization.

## 2.2 Natural frame for the L-curve during iterations

In the case of a poorly conditioned system, the reduction of the error can be obtained at the price of an explosion of the norm of the solution. This is well explained by Picard analysis: the phenomenon occurs when the eigenvalues of the operator decrease faster than the contribution of the right-hand side in the associated direction. It can also be visualized on a L-curve, in the positive quarter of a frame of the form ($\|\mathbf{r}_i\|$, $\|\mathbf{x}_i\|$): the curve starts in the bottom right corner (large error, small norm) with a fast decay of the error, and finishes in the top left corner (reduced error, large norm). As shown earlier, conjugate gradient provides natural norms to evaluate the error and the norm of the solution: $\|\mathbf{x}_i - \mathbf{x}_0\|_{\mathbf{A}}$ and $\|\mathbf{x}_i - \mathbf{x}_0\|_{\mathbf{M}}$. With this choice of norms, the curve is always oriented toward the upper-left corner: at each iteration, the norm of the error decreases and the norm of the solution increases.

## 2.3 A posteriori filtrering

Ritz elements offer a convenient way to filter the solution. Assuming $m$ iterations were conducted, we can process the basis $\mathbf{V}_m$ and the values $\Theta_m$. We can decompose the right-hand side on the Ritz basis $r_j^{(m)} = \mathbf{v}_j^{(m)T} \mathbf{r}_0$, and define:

$$ \text{for } i \le m, \quad \tilde{\mathbf{x}}_i^{(m)} = \mathbf{x}_0 + \sum_{j=1}^{i} \frac{r_j^{(m)}}{\theta_j^{(m)}} \mathbf{v}_j^{(m)}. \qquad (13) $$

We have:

$$ \|\tilde{\mathbf{x}}_i^{(m)} - \mathbf{x}\|_{\mathbf{A}}^2 = \|\tilde{\mathbf{x}}_i^{(m)} - \mathbf{x}_0\|_{\mathbf{A}}^2 - \sum_{j=1}^{i} \frac{(r_j^{(m)})^2}{\theta_j^{(m)}}, \qquad (14) $$

$$ \|\tilde{\mathbf{x}}_i^{(m)} - \mathbf{x}_0\|_{\mathbf{M}}^2 = \sum_{j=1}^{i} \frac{(r_j^{(m)})^2}{(\theta_j^{(m)})^2}, $$

and of course:

$$ \|\tilde{\mathbf{x}}_i^{(m)} - \tilde{\mathbf{x}}_{i-1}^{(m)}\|_{\mathbf{A}}^2 = \frac{(r_i^{(m)})^2}{\theta_i^{(m)}} \quad \text{and} \quad \|\tilde{\mathbf{x}}_i^{(m)} - \tilde{\mathbf{x}}_{i-1}^{(m)}\|_{\mathbf{M}}^2 = \frac{(r_i^{(m)})^2}{(\theta_i^{(m)})^2}. \qquad (15) $$

Since the $(\theta_i^{(m)})$ are sorted in decreasing order, we see that the error of $i \mapsto (\tilde{\mathbf{x}}_i^{(m)})$ tends to decrease slower than its norm tends to increase. The L-curve for $(\tilde{\mathbf{x}}_i^{(m)})_i$ is then convex and does not zig-zag.

The slope of the L-curve between the point $i-1$ and $i$ is $-(\theta_i^{(m)})^{-1}$. A possibility is to define the corner as the point which maximizes the variation of slope: $i = \arg\max_j((\theta_{j+1}^{(m)}) - (\theta_j^{(m)})^{-1})$.

Ritz' elements also make it possible to use Picard's theory and stop the construction of $\tilde{\mathbf{x}}_i^{(m)}$ when the contribution $j \mapsto r_j^{(m)}$ starts to decrease less fast than $j \mapsto \theta_j^{(m)}$. This criterion has the advantage to take into account the properties of the right-hand side.

# 3 Preconditioning by regularization

The core idea of this paper is that the preconditioner and the regularization should be the same operator. Thus, we use the matrix $\mathbf{M}$ for the Tikhonov regularization. As mentioned in the introduction, this idea makes sense as the same physical motivation underlies the choice of the regularization and that of the preconditioner. Moreover, many opportunities are opened by this choice. However, the hypothesis behind this idea is that there exists a cheap technique to apply the preconditioner (i.e. $\mathbf{M}^{-1}$).
As suggested by the notation, we investigate the effects of using the same operator for the regularization and the preconditioning, in particular when there exists a cheap technique to apply the preconditioner (i.e. $M^{-1}$). Conceptually, this idea makes sense as the same physical motivation underlies the choice of the regularization and that of the preconditioner. Moreover, many opportunities are opened by this choice.

If we assume that the system (16) was solved for a given $\lambda$ in $m$ iterations, then we can process the Ritz basis $V_m$. The strong point is that the properties of $V_m$ are independent of $\lambda$:

$$ 
\begin{aligned}
V_m^T M V_m &= I \\
V_m^T A_\lambda V_m &= \Theta_{\lambda,m} = \Theta_m + \lambda I.
\end{aligned}
\qquad (17) $$

*Remark 1.* $\lambda$ can be viewed as a shift in the generalized eigenvalues of $(A, M)$. Since $\lambda$ alters the initial residual and only a limited number of iterations is made, the content of $V_m$ is influenced by $\lambda$, without impairing the orthogonality properties.

Note that the initial residual takes the form:

$$ r_{\lambda,0} = b_{\lambda} - A_{\lambda}x_0 = \underbrace{(b_A - Ax_0)}_{r_{A,0}} + \lambda \underbrace{(b_M - Mx_0)}_{r_{M,0}} \qquad (18) $$

After $m$ iterations, we can define the Ritz' approximation:

$$ \tilde{\mathbf{x}}_{\lambda,i}^{(m)} = \mathbf{x}_0 + \sum_{j=1}^{i} \mathbf{v}_j^{(m)} \frac{(\mathbf{v}_j^{(m)})^T \mathbf{r}_{\lambda,0}}{\theta_j^{(m)} + \lambda} = \mathbf{x}_0 + \sum_{j=1}^{i} \mathbf{v}_j^{(m)} \frac{\mathbf{v}_j^{(m)T} \mathbf{r}_{A,0} + \lambda \mathbf{v}_j^{(m)T} \mathbf{r}_{M,0}}{\theta_j^{(m)} + \lambda} \qquad (19) $$

This approximation can be computed at zero cost, and its dependence in $\lambda$ is explicit: the L-curve of $\lambda \mapsto \tilde{\mathbf{x}}_{\lambda,i}^{(m)}$ can be plotted as a continuous function. It even permits to give sense to the limit solution when $\lambda \to 0$ even when $A$ was not invertible. It also gives an analytical formula for the search of the optimal choice of $(\lambda, i)$ realizing a good compromise between error and norm of the solution. Indeed, noting $r_{\lambda,j}^{(m)} = \mathbf{v}_j^{(m)T} \mathbf{r}_{\lambda,0}$, we have the properties:

$$ 
\begin{aligned}
\|\tilde{\mathbf{x}}_{\lambda,i}^{(m)} - \mathbf{x}_0\|_{\mathbf{M}}^2 &= \sum_{j=1}^{i} \left( \frac{r_{\lambda,j}^{(m)}}{\theta_j^{(m)} + \lambda} \right)^2 \\
\|\tilde{\mathbf{x}}_{\lambda,i}^{(m)} - \mathbf{x}_\lambda\|_{\mathbf{A}_\lambda}^2 &= \mathbf{r}_{\lambda,0}^T (\mathbf{A} + \lambda \mathbf{M})^{-1} \mathbf{r}_{\lambda,0} - \frac{\left( r_{\lambda,j}^{(m)} \right)^2}{\theta_j^{(m)} + \lambda}
\end{aligned}
\qquad (20) $$

Sadly, the measurement of the error in a $\lambda$-independent norm like $\|\tilde{\mathbf{x}}_{\lambda,i}^{(m)} - \mathbf{x}_\lambda\|_{\mathbf{A}}^2$ seems to be unavailable at the moment.

# 4 Assessments – Application to optical flow reconstruction

The recovery of the optical flow is briefly summarized in Appendix A. It is in fact a nonlinear minimization problem. A pyramidal approach is developed as often in image correlation, in order to provide a meaningful initialization [4]. For simplicity, we focus on the last nonlinear system to be solved, associated with the full image. Anyhow, the initialization of this system was impacted by the choice of the regularization.

We consider the solution to system (23, 24) with augmented preconditioned conjugate gradient, Algorithm 1. Due to the rectangular shape of the images, there exists an extremely cheap way to solve the preconditioner, which is a Laplace operator, using Fast Fourier transform or more precisely discrete cosine transform, see Appendix B.

The proposed test case is a holed composite plate in traction, with a 45° crack to be identified at the bottom of the hole. The speckle in the initial configuration is shown in Figure 1. To quantify the bad conditioning, the non-zero eigenvalues of $A$ are in the interval $[10^{-6}, 10^2]$.
We are interested in Tikhonov-regularized systems of the form:

$$ \underbrace{(\mathbf{A} + \lambda \mathbf{M})}_{\mathbf{A}_{\lambda}} \mathbf{x}_{\lambda} = \underbrace{\mathbf{b}_{\mathbf{A}} + \lambda \mathbf{b}_{\mathbf{M}}}_{\mathbf{b}_{\lambda}}. \qquad (16) $$

Note that we chose a $\lambda$-affine form for the right-hand side because it meets our practical needs, but the method applies to any separate form $(\mathbf{b}(\lambda) = \sum_a f_a(\lambda)\mathbf{b}_a)$.

If we assume that the system (16) was solved for a given $\lambda$ in $m$ iterations using the proposed M-preconditioned conjugate gradient algorithm, then we can process the Ritz basis $\mathbf{V}_m$. The strong point is that $\mathbf{V}_m$ separates the effects of the operator $\mathbf{A}$ and that of the regularization $\mathbf{M}$ independently of $\lambda$:

$$ \begin{aligned} \mathbf{V}_m^T \mathbf{M} \mathbf{V}_m &= \mathbf{I}_m, \\ \mathbf{V}_m^T \mathbf{A}_\lambda \mathbf{V}_m &= \mathbf{\Theta}_{\lambda,m} = \mathbf{\Theta}_m + \lambda \mathbf{I}_m. \end{aligned} \qquad (17) $$

*Remark 1.* $\lambda$ can be viewed as a shift in the generalized eigenvalues of $(\mathbf{A}, \mathbf{M})$. Since $\lambda$ alters the initial residual and only a limited number $m$ of iterations is made, the content of $\mathbf{V}_m$ depends on the choice of $\lambda$, but the orthogonality properties remain.

Note that the initial residual takes the separate form:

$$ \mathbf{r}_{\lambda,0} = \mathbf{b}_{\lambda} - \mathbf{A}_{\lambda} \mathbf{x}_{0} = \underbrace{(\mathbf{b}_{\mathbf{A}} - \mathbf{A} \mathbf{x}_{0})}_{\mathbf{r}_{\mathbf{A},0}} + \lambda \underbrace{(\mathbf{b}_{\mathbf{M}} - \mathbf{M} \mathbf{x}_{0})}_{\mathbf{r}_{\mathbf{M},0}}, \qquad (18) $$

and we can define the spectral contributions $r_{A,j}^{(m)} = \mathbf{v}_j^{(m)T} \mathbf{r}_{\mathbf{A},0}$ and $r_{M,j}^{(m)} = \mathbf{v}_j^{(m)T} \mathbf{r}_{\mathbf{M},0}$.

After $m$ iterations, we can define the Ritz' approximations for $i \le m$:

$$ \tilde{\mathbf{x}}_{\lambda,i}^{(m)} = \mathbf{x}_0 + \sum_{j=1}^{i} \mathbf{v}_j^{(m)} \frac{(\mathbf{v}_j^{(m)})^T \mathbf{r}_{\lambda,0}}{\theta_j^{(m)} + \lambda} = \mathbf{x}_0 + \sum_{j=1}^{i} \frac{r_{A,j}^{(m)} + \lambda r_{M,j}^{(m)}}{\theta_j^{(m)} + \lambda} \mathbf{v}_j^{(m)} \qquad (19) $$

These approximations can be computed at marginal cost, and there dependence in $\lambda$ is explicit: the L-curve of $\lambda \mapsto \tilde{\mathbf{x}}_{\lambda,i}^{(m)}$ is a rational fraction. It even permits to give sense to the limit solution when $\lambda \to 0$ even when $\mathbf{A}$ was not invertible. It also gives an analytical formula for the search of the optimal choice of $(\lambda, i)$ realizing a good compromise between error and norm of the solution. we have the properties:

$$ \begin{aligned} \|\tilde{\mathbf{x}}_{\lambda,i}^{(m)} - \mathbf{x}_0\|_{\mathbf{M}}^2 &= \sum_{j=1}^{i} \left( \frac{r_{A,j}^{(m)} + \lambda r_{M,j}^{(m)}}{\theta_j^{(m)} + \lambda} \right)^2 \\ \|\tilde{\mathbf{x}}_{\lambda,i}^{(m)} - \mathbf{x}_{\lambda}\|_{\mathbf{A}_{\lambda}}^2 &= \|\tilde{\mathbf{x}}_{\lambda,0}^{(m)} - \mathbf{x}_{\lambda}\|_{\mathbf{A}_{\lambda}}^2 - \sum_{j=1}^{i} \frac{\left(r_{A,j}^{(m)} + \lambda r_{M,j}^{(m)}\right)^2}{\theta_j^{(m)} + \lambda} \end{aligned} \qquad (20) $$

We write $\mathbf{V}_i^{(m)} = (\mathbf{v}_1^{(m)}, \dots, \mathbf{v}_i^{(m)})$ and $\mathbf{\Theta}_i^{(m)} = \mathrm{diag}(\theta_1^{(m)}, \dots, \theta_i^{(m)})$, so that $\tilde{\mathbf{x}}_{\lambda,i}^{(m)} = \mathbf{x}_0 + \mathbf{V}_i^{(m)}(\mathbf{\Theta}_i^{(m)} + \lambda \mathbf{I}_i)^{-1} \mathbf{V}_i^{(m)T} \mathbf{r}_{0,\lambda}$. We have an even more interesting measure of the error in non-regularized norm with respect to the non-regularized solution:

$$ \begin{aligned} \|\tilde{\mathbf{x}}_{\lambda,i}^{(m)} - \mathbf{x}\|_{\mathbf{A}}^2 - \|\mathbf{x}_0 - \mathbf{x}\|_{\mathbf{A}}^2 &= \|\mathbf{V}_i^{(m)}(\mathbf{\Theta}_i^{(m)} + \lambda \mathbf{I}_i)^{-1} \mathbf{V}_i^{(m)T} \mathbf{r}_{\lambda,0} + \mathbf{x}_0 - \mathbf{x}\|_{\mathbf{A}}^2 - \|\mathbf{x}_0 - \mathbf{x}\|_{\mathbf{A}}^2 \\ &= \|\mathbf{V}_i^{(m)}(\mathbf{\Theta}_i^{(m)} + \lambda \mathbf{I}_i)^{-1} \mathbf{V}_i^{(m)T} \mathbf{r}_{\lambda,0}\|_{\mathbf{A}}^2 \\ &\quad - 2\mathbf{r}_{\mathbf{A},0}^T \mathbf{V}_i^{(m)}(\mathbf{\Theta}_i^{(m)} + \lambda \mathbf{I}_i)^{-1} \mathbf{V}_i^{(m)T} \mathbf{r}_{\lambda,0} \\ &= \sum_{j=1}^{i} \frac{\left(r_{A,j}^{(m)} + \lambda r_{M,j}^{(m)}\right) \left(\theta_j^{(m)} \left(r_{A,j}^{(m)} + \lambda r_{M,j}^{(m)}\right) - 2r_{A,j}^{(m)}\right)}{(\theta_j^{(m)} + \lambda)} \end{aligned} \qquad (21) $$

Another feature is the possibility to analyze the system in terms of spectral content $(\theta_i^{(m)})$ vs regularization $\lambda$, and with regard to the decomposition of the right-hand side $(r_{A,j}^{(m)} + \lambda r_{M,j}^{(m)})$. This can be particularly well visualized in a Picard plot.
<figure><img src="image_3.png" /><figcaption>Figure 1: Speckle of the test specimen.</figcaption></figure>

## 4.1 Quality of the preconditioner

We first wish to verify that preconditioning by regularization actually leads to better enforcement of the regularity. In Table 1, we can qualitatively compare the classical Jacobi approach of preconditioning by the diagonal of the operator $M_{jac}^{-1} = \text{diag}(\mathbf{A}_{\lambda})^{-1}$ and the proposed preconditioning by regularization. The increased regularity is particularly visible for low weight $\lambda$ and low precision $\varepsilon$ of the linear solver.

Preconditioning by the regularization operator thus makes it possible to make meaningful computations with low weight in the regularization and to solve with less precision, hence with fewer iterations. Nevertheless, one has to mention that our preconditioner is computationally more expensive per iteration than the diagonal one.

An interesting scenario unfolds. The regularization preconditioner promotes low frequency corrections. Indeed, it is associated with a fully populated matrix (never actually computed) and the search directions have naturally large wavelength. As iterations progress, higher frequency modes emerge introducing more and more details and irregularity. On the contrary, the Jacobi preconditioner is diagonal, and it naturally encourages (independent) details, only iterations make it possible to reveal the structure between neighboring pixels.

## 4.2 Ritz filtering

We analyze the solving process for high ($\lambda = 1000$) and low ($\lambda = 10$) levels of regularization. We use the second stopping criterion of Equation 12 with $\varepsilon = 10^{-5}$, which corresponds to a rather high degree of convergence. The identified strain field are given in Figure 4.

We analyze the convergence in terms of compromise between the decrease of the error and the increase of the norm of the gradient of the solution which stems from the oscillations in the identified fields. Figure 2 presents two L-curves associated with high and low regularization. We use the natural CG-norms, please note that the position of the 0-abscissa is conventional because $\|\mathbf{x}_0 - \mathbf{x}\|_{\mathbf{A}}$ is unknown. The L-curves of the CG iterations (dotted lines) have similar shapes, like pieces of hyperbola. Due to the difference of magnitude, different scales had to be used: the error decreases four times less when the high regularization is used, and the norm of the solution remains 50 times smaller.

In order to better understand the convergence, we conduct a Ritz analysis. For $\lambda = 1000$, the convergence is attained in $m = 56$ iterations and as many Ritz vectors are computed. Table 2 presents a selection of these modes, sorted in decreasing order of Ritz value. The Ritz vectors resemble vibration modes with increasing number of anti-nodes. The first vectors are so regular that the hole is barely visible. The crack is only visible on the latest modes. This is bad (but logical) news because these are the most difficult modes to converge, thus they are probably bad approximations of actual eigenmodes, and the crucial mechanical information they carry is difficult to reuse.

We use formula (19) for the *a posteriori* filtering of the solution based on Ritz vectors. We present the L-curves in terms of modes included in the reconstruction. We show the curves in terms of full error and only taking into account the image error ($v_i^T r_{A,0}$) — they are almost overlaid on each other, a slight discrepancy only appears for high regularization. The shape
# 4 Assessment in the linear case: data completion problem

## 4.1 Laplace PDE with missing and redundant boundaries

We use the classical illustration of the ill-posedness for the inverse Laplace problem. We consider the rectangular domain $\Omega := [0, T] \times [0, H]$, where the following Cauchy problem holds:

$$ 
\begin{aligned}
\Delta u &= 0 && \text{in } \Omega, \\
u &= 0 && \text{on } y = 0 \text{ and } y = H, \\
\frac{\partial u}{\partial x} &= 0 && \text{on } x = 0, \\
u &= u_L := \sin k\pi \frac{y}{H} && \text{on } x = 0,
\end{aligned}
\quad (22) $$

where $k \in \mathbb{N}$ is the wavenumber of the signal. As can be observed there is no boundary condition on the right-hand side $\Gamma_R = \{(T, y), y \in (0, H)\}$ whereas there are both Dirichlet and Neumann conditions on the left-hand side. This is a model problem for the non-destructive control of structures in statics, it admits the following solution:

$$ u(x, y) = \sin k\pi \frac{y}{H} \cosh(k\pi \frac{x}{H}). \quad (23) $$

We observe that the solution on the right-hand side ($x = T$) explodes for thick domains ($T \nearrow$) or large wavenumbers ($k \nearrow$) as illustrated in Figure 1. This kind of problem is often qualified as “severely ill-posed” [3].

<figure><img src="image_2.png" /><figcaption>Figure 1: Solution to be identified by Steklov-Poincaré approach observing only the left-hand side (indiscernible oscillation), for $k = 3$ and $H = T = 1$.</figcaption></figure>

## 4.2 Steklov-Poincaré approach

This approach [4] relies on the reformulation of the Cauchy problem in two well-posed problems with one common unknown boundary value $u_R$ on the right-hand side and just one piece of information on the left-hand side.

$$ 
\begin{aligned}
\Delta u_D &= \Delta u_N = 0 && \text{in } \Omega, \\
u_D &= u_N = 0 && \text{on } y = 0 \text{ and } y = H, \\
u_D &= u_L \text{ and } \frac{\partial u_N}{\partial x} = 0 && \text{on } x = 0, \\
u_D &= u_N = u_R && \text{on } \Gamma_R.
\end{aligned}
\quad (24) $$

The index indicates whether Dirichlet or Neumann boundary conditions are considered on the left-hand side. Using classical variational theory, these problems have one solution in $H^1(\Omega)$. Using linearity, we can define the Steklov-Poincaré operators $H_{00}^{1/2}(\Gamma_R) \to H^{-1/2}(\Gamma_R)$:

$$ \frac{\partial u_D}{\partial x} = S_D(u_R) - b_D, \quad \frac{\partial u_N}{\partial x} = S_N(u_R), \quad \text{on } x = L. \quad (25) $$

$S_D$ and $S_N$ are linear operators and $b_D$ is the contribution of the non-zero Dirichlet condition $u_L$. Solving Cauchy problem (22) is then equivalent to finding $u_R$ such that:

$$ (S_D - S_N)u_R = b_D. \quad (26) $$
<table><thead><tr><th>λ</th><th>ε</th><th>Diagonal Prec.</th><th>Regularization Prec.</th></tr></thead><tbody><tr><td>low</td><td>low</td><td></td><td></td></tr><tr><td>low</td><td>high</td><td></td><td></td></tr><tr><td>high</td><td>low</td><td></td><td></td></tr><tr><td>high</td><td>high</td><td></td><td></td></tr></tbody></table>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Table 1: $\epsilon_{xx}$ strain field (range = mean value ± 3 st.dev.). Comparison of the effect of preconditioning by diagonal (simple approach) vs by regularization, for different weights $\lambda \in \{1, 1000\}$ and linear solver precision $\epsilon \in \{10^{-2}, 10^{-3}\}$.</figcaption></figure>

of the Ritz L-curves corresponds to most of the modes (the highest) only slightly decreasing the error and almost not changing the norm, only the last modes, which contain the crack information, are associated with significant decrease of the error (but of course at the cost of much increased solution norm). More or less, if a corner was to be selected it would correspond to just suppressing the contribution of the last mode.

In order to better understand this behavior, we first analyze the convergence of the Ritz values by comparing the spectrum obtained at the last iteration with the one obtained just one iteration before, like was done in [8] in the case of a well posed problem. It appears that the largest Ritz values were quite well approximated and only the lowest part of the spectrum evolves (remember the Ritz values correspond to the inverse of the slope of the segments in the L-curve). In other words, even though the last iterations seem not to modify the solution much (accumulation of the dots in the upper left part on the CG L-curves), they play an important role in terms of estimation of the lower part of the spectrum, without adding lots of small eigenvalues.

To support this analysis, we conducte a Picard's study on Figure 3 which shows the distribution of the Ritz values ($\theta_i^{(m)}$) as well as the decomposition of the right-hand side on the eigenspace ($v_i^{(m)T} r_{A,0}$) and $\lambda(v_i^{(m)T} r_{M,0})$. It is worth recalling that low and highly regularized systems have the same spectrum, except that it is more sampled for the low regularization which requires two times more iterations to converge. The Ritz values are slowly decreasing and only the last 10% really decay, the low regularization is not associated with an overpopulation of the lowest part of the spectrum. What stands out is the fact that the right-hand side contributes almost equally on all modes (at least it does not decrease for larger Ritz values). Picard's theory
The operator ($S_D - S_N$) is compact, meaning that its eigenvalues accumulate near zero. The problem can be discretized using finite element and solved iteratively. We use boldface for the discrete counterpart of the operators and fields introduced above. In [13] conjugate gradient with several preconditioners was tested. It was observed that trying to approximate the inverse of the linear operator was a bad idea whereas an efficient preconditioner was provided by the discrete $S_D$, which in fact corresponded using a Krylov solver to accelerate the KMF stationary iteration [22]. Note that $(u_R, S_D(u_R)) = \|u_D\|_{H^1(\Omega)}^2$ so that the KMF preconditioner is in fact a measure of the energy of the field.

## 4.3 Direct solvers

We consider the case $H = T = 1$ and $k = 3$, the domain is discretized using 40 × 40 square 4-node Lagrange elements. A Gaussian white noise is added to the signal with signal-to-noise ratio of 10 dB. The input signal is given in Figure 2a. Figure 2b presents the identified fields using either a truncated SVD where only singular values larger than $\epsilon_\sigma\sigma_1$ are kept, or a regularized operator ($S_D - S_N$) + $\lambda R$), $R$ being either identity $I$ or $S_D$. For the record, the five largest eigenvalues of ($S_D - S_N$) are {5.8 $10^{-4}$, 2.1 $10^{-6}$, 5.6 $10^{-9}$, 1.2 $10^{-11}$, 2.3 $10^{-14}$} while the spectrum of $S_D$ is contained in the interval [7.8 $10^{-3}$, 1.63].

<figure><img src="image_3.png" /></figure>

For all methods, the straight lines correspond to an identified value of the driving parameter which performs well, while the dotted lines correspond to too strong regularization, and the “+” lines correspond to insufficient regularization. We observe that the Tikhonov regularization by $S_D$ performs slightly better in the sense that the identification seems to be less sensitive to the value of the regularization parameter.

## 4.4 Preconditioned iterative solvers

We are now interested in solving the system with conjugate gradient and trying to find efficiently the optimal parameters for the best identification. We compare conjugate gradient without preconditioner $M = I$, with Jacobi preconditioner $M_J = \text{diag}(S_D - S_N)$ which is usually a good idea for well-posed problems, and the KMF preconditioner $M = S_D$. The stopping criteria is the one given in (12b) with $\epsilon = 10^{-9}$.

The first column of Figure 3 presents the successive approximations during CG (also compared with the reference and with the solution obtained with a truncated SVD with threshold $10^{-12}$). The first three rows correspond to non-regularized systems. It appears clearly that due to the quasi-singularity of the operator, the residual alone is not capable of providing a meaningful stopping criterion: for the first two rows the solvers does not stop at iteration 2 whereas the solution is quite close to the reference. The $S_D$ preconditioner appears not to behave significantly differently from the non-preconditioned case, which could be expected as its spectrum does not spread widely. Jacobi preconditioner (third row) behaves poorly, this is mostly due to the fact that the preconditioner is almost singular near the extremities. A fourth row was added which presents the same preconditioner with a slight Tikhonov regularization $\lambda S_D$ with $\lambda = 10^{-9}$, in that case iterations manage to reduce the error (note that the solver was stopped after 20 iterations).

The second and third columns present the iteration L-curves using either the Euclidean frame ($\|r\|_2^2$, $\|x_i - x_0\|_2^2$) or the natural frame ($\|x_i - x\|_\infty$, $\|x_i - x_0\|_\infty$). As could be feared, the Euclidean frame gives hard to exploit curves where it is impossible to define a corner. The natural frame provides much nicer convex curves, and at least in the first two rows, one can clearly see that the last iteration barely reduces the error while strongly increasing the norm of the solution, which suggests selecting the solution at the second to last iteration (which agrees with the “eye norm” on the first column).
<figure><img src="image_12.png" /><figcaption>Table 2: $\epsilon_{xx}$ strain field for 8 Ritz vectors out of the 56 computed ($\epsilon = 10^{-5}$, $\lambda = 1000$).</figcaption></figure>

thus suggests that we should stop the reconstruction when the Ritz values start to decay. This is not possible in our case since the crack is mostly represented in this part of the spectrum.

By the way, Figure 3 permits to compare the smallest Ritz value $\theta_{\min}$ with the regularization parameter $\lambda$. The case that we called “low regularization” corresponds to $\lambda$ being negligible with respect to the small Ritz value $\theta_m^{(m)}$, and thus only marginally modifying the active Ritz spectrum. On the contrary, the high regularization corresponds to a $\lambda > \theta_m^{(m)}$ which means that the lower part of the spectrum of $A_{\lambda}$ is flattened relative to that of $A$.

<figure><img src="image_13.png" /><figcaption>Figure 2: L-curves (same $\epsilon = 10^{-5}$) compared with Ritz post-treatment, for different regularization intensity $\lambda$. Note that different scales are used on the plots.</figcaption></figure>

## 4.3 Subspace recycling

Even though it appears that Ritz filtering is difficult to apply to the studied system, we can still benefit from Ritz vectors to accelerate the solution. As a sequence of linear systems with identical matrix has to be solved, it is natural to augment the system with the previously generated Ritz vectors by concatenating $C \leftarrow (C V_m)$. Indeed, augmentation comes with optimized block operations that make augmenting by one vector much cheaper than one iteration.

Moreover, Ritz vectors possess two advantages. Firstly, the product $AV_m$ which is required during augmentation can be obtained at low computational cost using the formula: $A\hat{z}_{j+1} = (-1)^{j+1}(q_{j+1} - \beta_j q_j)/\sqrt{\gamma_{j+1}}$ and $AV_m = A\hat{Z}_m\Xi_m$. Secondly, using normalization $V_m \leftarrow V_m\Theta_m^{-1/2}$, we have $V_m^TAV_m = I$.

Table 3 illustrates the performance of recycling for the nonlinear system to be solved on the
<table><caption>Figure 3: Data completion problem: quality of solution using CG with various preconditioner – effect of slight Tikhonov regularization on the worst preconditioner.</caption><thead><tr><th></th><th>prec.</th><th>Iterations</th><th>Euclidean L-curve</th><th>Natural L-curve</th></tr></thead><tbody><tr><td>Jacobi λ = 10<sup>-9</sup></td><td>KMF λ = 0</td><td></td><td></td><td></td></tr><tr><td>Jacobi λ = 0</td><td>KMF λ = 0</td><td></td><td></td><td></td></tr><tr><td>Jacobi λ = 0</td><td>KMF λ = 0</td><td></td><td></td><td></td></tr></tbody></table>
<figure><img src="image_14.png" /></figure>

$$ (a) \lambda = 1000 \approx \min(\theta_j) $$

<figure><img src="image_15.png" /></figure>

$$ (b) \lambda = 10 \ll (\theta_j) $$

<figure><img src="image_16.png" /><figcaption>Figure 3: Spectral analysis of the system for $\epsilon = 10^{-5}$ and different regularization intensity $\lambda$. A 5-width median filter was used to smooth out the contribution curves.</figcaption></figure>

$$ (a) \lambda = 1000 $$

<figure><img src="image_17.png" /></figure>

$$ (b) \lambda = 10 $$

<caption>Figure 4: Identified $\epsilon_{xx}$ field, for different regularization intensity $\lambda$, with $\epsilon = 10^{-5}$.</caption>

full image at the end of identification. The first linear system, only augmented by the kernel of the preconditioner is solved in 77 iterations. Then a certain portion of the Ritz vectors is used to augment the next 8 linear systems (same matrix, different right-hand sides). As the augmentation results in an excellent initialization, we use a criterion in terms of absolute value of $\|\mathbf{r}_j\|_M$ to halt the iterations because other comparison as given in Equation 12 might use an unfair reference. We measure the performance in terms of gain in iterations, and in computational time (measures are conducted on a upper mid-range laptop with Nvidia RTXA2000 graphic card). The gain in terms of iterations in moderate, with best obtained for small augmentation space (at most 1.3 iterations per augmentation vector, for 10 vectors). In terms of time, the optimal is obtained for augmentation space of 80%-90% of available vectors, with a global CPU time divided by almost 2 (this time includes all the extra cost associated with computing and using Ritz vectors). This size of subspace agrees with what we observed on the stability of the largest Ritz values in the L-curves plots.

## 4.4 Tuning of $\lambda$

It is often hard to automatize the selection of the regularization intensity $\lambda$. Picard's plots like in Figure 3 permit to put $\lambda$ in relation with the spectrum of the preconditioned operator and thus to understand the effect of the regularization in terms of flattened spectrum. Still, the final judge is often the expert's impression of a strain map, and it is convenient to compute maps associated with several ($\lambda_k$) at low cost.

Formula (19) makes it possible, after the solution of one linear system for a given $\lambda_0$, to post-
## 4.5 Opportunities when preconditioning with Tikhonov operator

We now consider the KMF operator $S_D$ used as both the regularization and the preconditioner: $(A+\lambda M)x = b$. After the $m$ CG iterations needed to reach convergence, we postprocess the Ritz values $\Theta^{(m)}$, and the Ritz vectors $V^{(m)}$.

Figure 4 presents the approximation of the $\lambda$-L-curve given by the postprocessed Ritz' elements (21) from one iterative computation with $\epsilon = 10^{-9}$ and $\lambda = 10^{-9}$. It is compared to the L-curve obtained from direct solves with various $\lambda \in [10^{-6}, 10^{-12}]$, using the fact that in this academic case the reference $x$ is known. Note that because for the iterative solver the actual error is known up to the additive factor $\|x - x_0\|_A^2$, the curves are in fact translated so that their vertical asymptote is aligned with the y-axis.

We see that the L-curve postprocessed at negligible cost after the iterative solution provides an excellent estimation of the actual L-curve which would not be computable in real cases.

<figure><img src="image_4.png" /><figcaption>Figure 4: L-curve obtained with direct solves vs L-curve estimated by Ritz' approximation (21) after one iterative resolution with $\lambda = 10^{-9}$ materialized by the dot ($\epsilon = 10^{-9}$). The vertical asymptotes are aligned with the y-axis.</figcaption></figure>

Figure reffig:picardstek presents the Picard plot obtained after the iterative resolution with $\lambda = 10^{-9}$ and $\epsilon = 10^{-12}$ (for this simple case, higher precision was required in order to have more points to plot). First, we can observe the spectrum of $A$ estimated by the Ritz values. Thanks to (17), the regularization is a horizontal line, and we can directly see how the regularization impacts the spectrum: in this case, the two smallest eigenvalues are rectified.

A second interest of Picard plot is to compare the spectrum with the contribution of the right-hand side $|r_j^{(m)}|$. Clearly, beside the first two components, the right-hand side decreases slower than the Ritz' values. This indicates that the system is unstable as the effect of the contributions to the solution gets increasingly larger.

These two pieces of information observable in Picard plots enable informed choices of $\lambda$ and of the truncation index $i$ in the Ritz reconstruction $\tilde{x}_{\lambda,i}^{(m)}$.

# 5 Augmented preconditioned conjugate gradient for sequences of regularized systems

We now consider of nonlinear problems whose solution can be obtained by solving a sequence of linear systems with constant left-hand side and varying right-hand side:

$$ \text{For given } \lambda, \text{ at outer iteration } k \left\{ \begin{array}{l} \text{Compute right-hand side } \mathbf{b}_{\lambda}^{k}(\mathbf{x}_{\lambda}^{k}) \\ \text{Solve } \mathbf{A}_{\lambda}\boldsymbol{\delta}_{\lambda}^{k} = \mathbf{b}_{\lambda}^{k} \\ \text{Update } \mathbf{x}_{\lambda}^{k+1} = \mathbf{x}_{\lambda}^{k} + \boldsymbol{\delta}_{\lambda}^{k} \end{array} \right. \qquad (27) $$

To fully benefit the repetitive nature of the solves, we introduce augmentation in the conjugate gradient. Augmentation is a technique where the search space is split in two subspaces: the augmentation space, defined
<table><caption>Table 3: Performance of recycling</caption><tr><td>Aug.</td><td>0</td><td>10</td><td>20</td><td>30</td><td>40</td><td>50</td><td>60</td><td>70</td><td>max (77)</td></tr><tr><td>Iter.</td><td>77</td><td>64</td><td>57</td><td>49</td><td>44</td><td>41</td><td>40</td><td>40</td><td>38</td></tr><tr><td>Time (s)</td><td>11.7</td><td>10.1</td><td>8.3</td><td>7.1</td><td>6.6</td><td>6.3</td><td>6.2</td><td>6.3</td><td>6.6</td></tr></table>

<figure><img src="image_18.png" /><figcaption>Figure 5: Costless postprocessing for different ($\lambda_i$): top, initial computation with $\lambda_0 \in \{10, 100, 1000\}$; bottom, direct solution with $\lambda = 1$ and solutions deduced with $\lambda_1 = 1$. $\epsilon_{xx}$ strain field.</figcaption></figure>

process the solution for any $\lambda_k$ at the simple cost of computing the associated right-hand side $b_{A\lambda_k}$ (which depends on the history of the nonlinear solution for $\lambda_k$), and basic linear algebra operations.

Figure 5 presents the solution deduced for $\lambda_1 = 1$ from initial computations with different $\lambda_0$ (in $\{10, 100, 1000\}$) and $\epsilon = 10^{-4}$. Again, a median filter was applied after the Ritz reconstruction. It seems that the Ritz vectors make it possible to postprocess a reasonable solution with a $\lambda_1 \simeq \lambda_0/1000$. The reconstructed strain field appears to be much less smooth than the original computation (with $\lambda_0$) while less noisy than the direct low-regularization computation with $\lambda = 1$. If the deduced solution is not fully satisfying, it can still be used as an excellent initialization for a regular computation.

# 5 Conclusion

In this paper, we have studied how preconditioning and Tikhonov regularization could be efficiently combined in an augmented preconditioned conjugate gradient. We have shown that this association makes sense from a physical point of view and it made it possible to combine filtering, recycling of subspaces, and postprocessing of all regularized solutions at zero cost. This gives a favorable framework to apply criteria like the L-curve or Picard's analysis.

The solver was applied to a problem of optical flow reconstruction which introduced the extra difficulty of nonlinearity and the fact that the most important information was buried in the lower part of the spectrum. Satisfying results were obtained on actual measurements from digital image correlation of a mechanical test. Postprocessed solutions at zero-cost were still relevant for regularization weight $\lambda$ divided by up to 1000.

An obvious next step for this work is to consider inexact preconditioners, that is to say when the $M^{-1}$ matrix in the preconditioning step of the algorithm is only an approximation of the inverse of the regularization matrix $M$ in the operator. This would make the method applicable on a much broader class of problems.
<figure><img src="image_5.png" /><figcaption>Figure 5: Picard plot processed after the iterative resolution with $\lambda = 10^{-9}$ and $\varepsilon = 10^{-12}$.</figcaption></figure>

as the range of the given $n \times n_C$ full-rank matrix $C$, and its $A$-orthogonal complementary subspace. The part of the solution in the augmentation space is obtained at the initialization whereas the complement is searched for iteratively while maintaining the residual orthogonal to the augmentation.

We introduce the augmented Krylov subspace $\mathcal{K}_i(\mathbf{M}^{-1}\mathbf{A}, \mathbf{C}, \mathbf{M}^{-1}\mathbf{r}_0)$ [10]:

$$ \mathcal{K}_i(\mathbf{M}^{-1}\mathbf{A}, \mathbf{C}, \mathbf{M}^{-1}\mathbf{r}_0) = \operatorname{span} \left( \mathbf{M}^{-1}\mathbf{r}_0, \dots, (\mathbf{M}^{-1}\mathbf{A})^{(i-1)}\mathbf{M}^{-1}\mathbf{r}_0 \right) \oplus \operatorname{Range}(\mathbf{C}) \quad (28) $$

Given an arbitrary initialization $\mathbf{x}_{00}$ and associated residual $\mathbf{r}_{00} = \mathbf{b} - \mathbf{A}\mathbf{x}_{00}$, the $i_{th}$ iteration can be defined as:

$$ \left\{ \begin{array}{ll} \text{find} & \mathbf{x}_i \in \mathbf{x}_{00} + \mathcal{K}_i(\mathbf{M}^{-1}\mathbf{A}, \mathbf{C}, \mathbf{M}^{-1}\mathbf{r}_{00}) \\ \text{such that} & \mathbf{r}_i \perp \mathcal{K}_i(\mathbf{M}^{-1}\mathbf{A}, \mathbf{C}, \mathbf{M}^{-1}\mathbf{r}_{00}) \end{array} \right. \quad (29) $$

This iteration is achieved by Algorithm 2 where the augmentation is managed by the correction of the initialization (in order to obtain $\mathbf{x}_0$) and the projector $\mathbf{P}$ on $\ker(\mathbf{C}^T\mathbf{A})$, which together ensure that the residual remains orthogonal to $\operatorname{Range}(\mathbf{C})$ [12].

<caption>**Algorithm 2 Augmented Preconditioned Conjugate Gradient**</caption>

x₀₀ and C given
P = I - C(C^TAC)^{-1}C^T A
x₀ = P x₀₀ + C(C^TAC)^{-1}C^T b = x₀₀ + C(C^TAC)^{-1}C^T r₀₀
r₀ = b - A x₀₀ = P^T r₀₀
z₀ = P M^{-1} r₀₀, w₀ = z₀
γ₀ = (z₀^T r₀)
for i = 0, 1, ..., m (convergence) do
    qᵢ = A wᵢ
    δᵢ = (wᵢ^T qᵢ), αᵢ = δᵢ^{-1} γᵢ
    xᵢ₊₁ = xᵢ + wᵢ αᵢ
    rᵢ₊₁ = rᵢ - qᵢ αᵢ
    zᵢ₊₁ = P M^{-1} rᵢ₊₁
    γᵢ₊₁ = (zᵢ₊₁^T rᵢ₊₁)
    βᵢ = γᵢ^{-1} γᵢ₊₁
    wᵢ₊₁ = zᵢ₊₁ + wᵢ βᵢ
end for

## 5.1 Roles of augmentation

One particular use of augmentation is to handle symmetric positive semi-definite preconditioner $\mathbf{M}$: by ensuring $\ker(\mathbf{M}) \subset \operatorname{Range}(\mathbf{C})$ we make sure that the (pseudo-)inverse of $\mathbf{M}$ is well-defined for the residuals it is applied to. Let $\mathbf{C}_0$ be a basis of $\ker(\mathbf{M})$.
References

[1] Owe Axelsson and Igor Kaporin. Error norm estimation and stopping criteria in preconditioned conjugate gradient iterations. *Numerical Linear Algebra with Applications*, 8(4):265-286, 2001.

[2] Owe Axelsson and Gunhild Lindskog. On the rate of convergence of the preconditioned conjugate gradient method. *Numerische Mathematik*, 48:499-523, 1986.

[3] Gilles Besnard, François Hild, and Stéphane Roux. Finite-element displacement fields analysis from digital images: application to portevin-le châtelier bands. *Experimental Mechanics*, 46:789-804, 2006.

[4] Ahmed Chabib, Jean-François Witz, Pierre Gosselet, and Vincent Magnier. The impact of metrics in mechanical imaging. preprint, 2023.

[5] Andrew Chapman and Youssef Saad. Deflated and augmented Krylov subspace techniques. *Numerical Linear Algebra with Applications*, 4(1):43-66, 1997.

[6] Zdeněk Dostál. Conjugate gradient method with preconditioning by projector. *International Journal of Computer Mathematics*, 23:315-323, 1988.

[7] Renaud Ferrier, Mohamed L. Kadri, and Pierre Gosselet. The Steklov-Poincaré technique for data completion: Preconditioning and filtering. *International Journal for Numerical Methods in Engineering*, 116(4):270-286, 2018.

[8] Pierre Gosselet, Christian Rey, and Julien Pebrel. Total and selective reuse of Krylov subspaces for the resolution of sequences of nonlinear structural problems. *International Journal for Numerical Methods in Engineering*, 94(1):60-83, 2013.

[9] Per Christian Hansen. The truncated SVD as a method for regularization. *BIT Numerical Mathematics*, 27(4):534-553, 1987.

[10] Per Christian Hansen. The discrete Picard condition for discrete ill-posed problems. *BIT Numerical Mathematics*, 30(4):658-672, 1990.

[11] Per Christian Hansen. Analysis of discrete ill-posed problems by means of the L-curve. *SIAM review*, 34(4):561-580, 1992.

[12] Zhongxiao Jia and G.W. Stewart. On the convergence of the Ritz values, Ritz vectors and refined Ritz vectors. Technical Report 3896, Institute of Advanced Computer Studies, Department of Computer Science, University of Maryland at College Park, 1999.

[13] Mohamed Larbi Kadri, Jalel Ben Abdallah, and Thouraya Nouri Baranger. Identification of internal cracks in a three-dimensional solid body via Steklov-Poincaré approaches. *Comptes Rendus Mécanique*, 339(10):674-681, 2011.

[14] Louis Kovalevsky and Pierre Gosselet. A quasi-optimal coarse problem and an augmented Krylov solver for the Variational Theory of Complex Rays. *International Journal for Numerical Methods in Engineering*, 2015.

[15] Vladimir Alekseevich Morozov. The error principle in the solution of operational equations by the regularization method. Zhurnal Vychislitel'noi Matematiki i Matematicheskoi Fiziki, 8(2):295-309, 1968.

[16] Youssef Saad. *Iterative methods for sparse linear systems*. SIAM, 2003.

[17] Youssef Saad. *Numerical Methods for Large Eigenvalue Problems*, volume 66 of *Classics in Applied Mathematics*. SIAM, Philadelphia, USA, revised edition, 2011.

[18] Andreǐ Nikolaevich Tikhonov and Vasiliǐ Iakovlevich Arsenin. *Solutions of ill-posed problems*. Vh Winston, 1977.
An additional use of augmentation is to reuse numerical information generated during one solve to accelerate the following. Indeed, augmentation comes with optimized block operations that make augmenting by one vector much cheaper than one iteration.

In our case, the first solve (in $m$ iterations) will serve to postprocess the Ritz basis $V_m$ and the next solves will be augmented by the concatenation $C = [C_0, V_{\tilde{m}}]$ where $\tilde{m} \le m$ indicates that only a selection of the first Ritz vectors is used. Indeed, storing too many vectors can result in excessive memory usage, moreover, the first Ritz vectors are often better converged and improve convergence further than the last [14].

Moreover, Ritz vectors possess two advantages. Firstly, the product $AV_m$ which is required during augmentation can be obtained at low computational cost using the formula:

$$AV_m = A\hat{Z}_m\Xi_m, \qquad (30) \\ \text{and } A\hat{Z}_{j+1} = (-1)^{j+1}(q_{j+1} - \beta_j q_j)/\sqrt{\gamma_{j+1}}.$$

Secondly, using normalization:

$$V_m \leftarrow V_m \Theta_m^{-1/2} \text{ leads to } V_m^T A V_m = I. \qquad (31)$$

## 5.2 Postprocessing of multiple solutions for various regularization

The nonlinearity makes it impossible to compute a full analytical L-curve as with formula (21). Nevertheless, we can choose a family of $P$ regularization coefficients $(\lambda_p)_{p \le P}$ for which we can obtain a good approximation at low cost.

$$\text{At outer iteration } k \left\{ \begin{array}{l} \text{Compute right-hand sides } \mathbf{b}_{\lambda_p}^k(\mathbf{x}_{\lambda_p}^k), \text{ for all } 0 \le p \le P \\ \text{Solve } A_{\lambda_0} \delta_{\lambda_0}^k = \mathbf{b}_{\lambda_0}^k \text{ extract Ritz basis } V_{m_k}, \\ \text{Estimate } \tilde{\delta}_{\lambda_p}^{k,(m_k)} \text{ for all } 0 < p \le P \text{ using (19)}, \\ \text{Update } \mathbf{x}_{\lambda_p}^{k+1} = \mathbf{x}_{\lambda_p}^k + \tilde{\delta}_{\lambda_p}^{k,(m_k)}, \text{ for all } 0 \le p \le P. \end{array} \right. \qquad (32)$$

In words, at each outer iteration, only the system associated with $\lambda_0$ is solved with APCG whereas the solutions for the remaining coefficients $(\lambda_p)_{0<p}$ are obtained thanks to Ritz' approximation.

# 6 Assessment in the nonlinear case: recovery of the optical flow

## 6.1 Optical flow in a nutshell

The optical flow is a digital image correlation technique which aims at estimating the displacement field between two images at the scale of the pixel. Contrarily to very popular approaches in solid mechanics inspired by the Finite Element Method [6], it does not rely on a mesh and on shape functions to approximate the displacement field. Given a sequence of two images $(I_1, I_2)$, viewed as $N \times M$ arrays of gray level pixels (in the discrete segment $\{0, 1, \dots, G_{max}\}$), it directly aims at finding the transformation $\phi = (\phi_x, \phi_y)$ such that $I_1 - I_2 \circ \phi = 0$. Note that we use interpolation between pixels so that the images can be defined on the rectangle $[0, N] \times [0, M] \subset \mathbb{R}^2$ with values in the continuous segment $[0, G_{max}] \subset \mathbb{R}$ and the displacement $u := (\phi - \mathcal{I})$ can take non-integer values ($\mathcal{I}$ is the identity operator). It is even common to obtain precision below one tenth of a pixel. In order to gain flexibility, and adapt to unavoidable noisy measurements which make the zero unachievable, the problem is better rephrased in terms of the minimization of the “image energy” $E_I$:

$$E_I^2 = \frac{1}{2} \|I_1 - I_2 \circ \phi\|^2, \quad \text{where } \|I\|^2 = \sum_{\substack{0 \le i < N \\ 0 \le j < M}} I(i, j)^2. \qquad (33)$$

Even under that form the problem is not well-posed, would it only be because there are two times more unknowns than equations. A solution to recover a well-posed problem is to enforce regularity to the displacement field. A penalty term related to the gradient is then introduced:

$$E^2 = E_I^2 + \frac{\lambda}{2} \|\nabla u\|^2, \qquad (34)$$

where we kept the Euclidean norm notation for $\|\nabla u\|^2 := \|\partial_x u_x\|^2 + \|\partial_y u_y\|^2 + \|\partial_x u_y\|^2 + \|\partial_y u_x\|^2$. $\lambda$ is a weight that needs to be tuned in order to balance the contributions of the image energy and of the regularization. Note that the quadratic minimization framework employed here is not necessarily the most relevant in terms of quality of the identified fields [30, 9] and that iterative methods were also developed for more sophisticated metrics [26].
A Basic notions about the optical flow

The optical flow is a digital image correlation technique which aims at estimating the displacement field between two images at the scale of the pixel. Contrarily to very popular approaches in solid mechanics inspired by the Finite Element Method [3], it does not rely on a mesh and on shape functions to approximate the displacement field. Given a sequence of two images ($I_1, I_2$), viewed as $N \times M$ arrays of gray level pixels (in the discrete segment $\{0, 1, \dots, G_{max}\}$), it directly aims at finding the transformation $\phi = (\phi_x, \phi_y)$ such that $I_1 - I_2 \circ \phi = 0$. Note that we use interpolation between pixels so that the images can be defined on the rectangle $[0, N] \times [0, M] \subset \mathbb{R}^2$ with values in the continuous segment $[0, G_{max}] \subset \mathbb{R}$ and the displacement $u := (\phi - I_d)$ can take non-integer values ($I_d$ is the identity operator). It is even common to obtain precision below one tenth of a pixel. In order to gain flexibility, and adapt to unavoidable noisy measurements which make the zero unachievable, the problem is better rephrased in terms of the minimization of the “image energy” $E_I$:

$$E_I^2 = \frac{1}{2} \|I_1 - I_2 \circ \phi\|^2, \quad \text{where } \|I\|^2 = \sum_{\substack{0 \le i < N \\ 0 \le j < M}} I(i, j)^2. \qquad (21)$$

Even under that form the problem is not well-posed, would it only be because there are two times more unknowns than equations. A solution to recover a well-posed problem is to enforce regularity to the displacement field. A penalty term related to the gradient is then introduced:

$$E^2 = E_I^2 + \frac{\lambda}{2} \|\nabla u\|^2, \qquad (22)$$

where we kept the Euclidean norm notation for $\|\nabla u\|^2 := \|\partial_x u_x\|^2 + \|\partial_y u_y\|^2 + \|\partial_x u_y\|^2 + \|\partial_y u_x\|^2$. $\lambda$ is a weight that needs to be tuned in order to balance the contributions of the image energy and of the regularization.

In general, a modified Newton approach is used to minimize the energy. Starting from a guess $u$, the update $u + du$ is computed by solving the system:

$$(\mathbf{A} + \lambda \mathbf{M})\mathbf{x} = \mathbf{b}_{\mathbf{A}} + \lambda \mathbf{b}_{\mathbf{M}}, \qquad (23)$$

with

$$\mathbf{A} = \begin{pmatrix} \mathbf{J}_x & \\ & \mathbf{J}_y \end{pmatrix} \begin{pmatrix} \mathbf{I} & \mathbf{I} \\ \mathbf{I} & \mathbf{I} \end{pmatrix} \begin{pmatrix} \mathbf{J}_x & \\ & \mathbf{J}_y \end{pmatrix}, \quad \mathbf{M} = \begin{pmatrix} \mathbf{\Delta} \\ \mathbf{\Delta} \end{pmatrix} \\ \mathbf{x} = \begin{pmatrix} \mathrm{vec}(du_x) \\ \mathrm{vec}(du_y) \end{pmatrix}, \quad \mathbf{b}_{\mathbf{A}} = \begin{pmatrix} \mathrm{vec}((I_1 - I_2 \circ \phi)J_x) \\ \mathrm{vec}((I_1 - I_2 \circ \phi)J_y) \end{pmatrix}, \quad \mathbf{b}_{\mathbf{M}} = \begin{pmatrix} \mathrm{vec}(\Delta u_x) \\ \mathrm{vec}(\Delta u_y) \end{pmatrix}. \qquad (24)$$

The vec operator converts images to vectors ($N \times M$ array to $NM$ vector). For $z \in \{x, y\}$, $J_z$ is the $z$ component of the gradient of $I_1$, $\Delta u_z$ is the (scalar) Laplace operator applied to $u_z$. $\mathbf{J}_z$ is the $NM$ diagonal operator containing the values of the gradient $J_z$. $\mathbf{I}$ and $\mathbf{\Delta}$ are respectively the $NM$ identity matrix and the $NM$ matrix version of Laplace operator (with Neumann boundary conditions). All the operators are in fact obtained by discrete difference on the image. Note that the gradient of $I_1$ is used to approximate the current Jacobian. As commonly done in image treatment, a median filter is applied to all the computed increments in order to remove outliers caused by the imperfect speckle.

It is extremely simple to work with $\mathbf{A}$ and $\mathbf{M}$ without assembling them, one only needs to compute and store the two $N \times M$ images $(\mathbf{J}_x, \mathbf{J}_y)$ and use Hadamard product and Laplace function when computing matrix-vector multiplication.

The system is of dimension $2MN$. As said earlier, $\mathbf{A}$ is strongly deficient since its rank is at most $MN$, a first part of its kernel has the following basis:

$$\mathrm{span} \begin{pmatrix} \mathbf{J}_y \\ -\mathbf{J}_x \end{pmatrix} \subset \ker(\mathbf{A}). \qquad (25)$$
A modified Gauss-Newton approach is used to minimize the energy [27] which we combine with a pyramidal approach in order to provide meaningful initializations (a sequence of reduced systems is defined, and starting from the coarsest the solution obtained at one level is extrapolated on the next level to define a sound initialization). For a given level of the pyramid, starting from a guess $u$, the update $u + du$ is computed by solving the system:

$$ (\mathbf{A} + \lambda \mathbf{M})\mathbf{x} = \mathbf{b}_{\mathbf{A}} + \lambda \mathbf{b}_{\mathbf{M}}, \qquad (35) $$

with

$$ \mathbf{A} = \begin{pmatrix} \mathbf{J}_x & & & \\ & \mathbf{J}_y & & \\ & & \mathbf{I} & \\ & & & \mathbf{I} \end{pmatrix} \begin{pmatrix} \mathbf{J}_x & & & \\ & \mathbf{J}_y & & \\ & & \mathbf{I} & \\ & & & \mathbf{I} \end{pmatrix}, \quad \mathbf{M} = \begin{pmatrix} \mathbf{\Delta} \\ & \mathbf{\Delta} \end{pmatrix} \\ \mathbf{x} = \begin{pmatrix} \mathrm{vec}(du_x) \\ \mathrm{vec}(du_y) \end{pmatrix}, \quad \mathbf{b}_{\mathbf{A}} = \begin{pmatrix} \mathrm{vec}((I_1 - I_2 \circ \phi)J_x) \\ \mathrm{vec}((I_1 - I_2 \circ \phi)J_y) \end{pmatrix}, \quad \mathbf{b}_{\mathbf{M}} = \begin{pmatrix} \mathrm{vec}(\Delta u_x) \\ \mathrm{vec}(\Delta u_y) \end{pmatrix}. \qquad (36) $$

The vec operator converts images to vectors ($N \times M$ array to $NM$ vector). For $z \in \{x, y\}$, $J_z$ is the $z$ component of the gradient of $I_1$, $\Delta u_z$ is the (scalar) Laplace operator applied to $u_z$. $\mathbf{J}_z$ is the $NM$ diagonal operator containing the values of the gradient $J_z$. $\mathbf{I}$ and $\mathbf{\Delta}$ are respectively the $NM$ identity matrix and the $NM$ matrix version of Laplace operator (with Neumann boundary conditions). All the operators are in fact obtained by discrete difference on the image. Note that the gradient of $I_1$ is used to approximate the current Jacobian. As commonly done in image treatment, a median filter is applied to all the computed increments in order to remove outliers caused by the imperfect speckle.

The proposed test case is a holed composite plate in traction, with a 45° crack to be identified at the bottom of the hole. The speckle in the initial configuration is shown in Figure 6. For the illustrations, we present the strain field obtained by deriving the computed displacement. Strain is indeed the mechanical quantity of interest to identify the crack. Only the $xx$ component is given as other components do not provide any further information. To quantify the bad conditioning, the non-zero eigenvalues of $\mathbf{A}$ are in the interval $[10^{-6}, 10^2]$.

<figure><img src="image_6.png" /><figcaption>Figure 6: Speckle of the test specimen.</figcaption></figure>

The test case as well as the method (with a simpler linear solver) are presented in [8]. For simplicity, we focus on the last nonlinear system to be solved, associated with the full image. Anyhow, the initialization of this system, resulting from the pyramidal approach, was impacted by the choice of the regularization.

## 6.2 Matrix free implementation

We consider the solution to system (35, 36) with preconditioned conjugate gradient.

Due to the rectangular shape of the images, there exists an extremely cheap way to solve the regularization matrix, which is a Laplace operator, using Fast Fourier transform or more precisely discrete cosine transform, see Appendix A.

It is extremely simple to work with $\mathbf{A}$ and $\mathbf{M}$ without assembling them, one only needs to compute and store the two $N \times M$ images $(\mathbf{J}_x, \mathbf{J}_y)$ and use Hadamard product and Laplace function when computing matrix-vector multiplication.

The system is of dimension $2MN$. As said earlier, $\mathbf{A}$ is strongly deficient since its rank is at most $MN$, a first part of its kernel has the following basis:

$$ \mathrm{span} \left( \begin{pmatrix} \mathbf{J}_y \\ -\mathbf{J}_x \end{pmatrix} \right) \subset \ker(\mathbf{A}). \qquad (37) $$

The rest of the spectrum is easy to compute since:

$$ \mathbf{A} \begin{pmatrix} \mathbf{J}_x \\ \mathbf{J}_y \end{pmatrix} = \begin{pmatrix} \mathbf{J}_x \\ \mathbf{J}_y \end{pmatrix} (\mathbf{J}_x^2 + \mathbf{J}_y^2). \qquad (38) $$
The rest of the spectrum is easy to compute since:

$$A \begin{pmatrix} J_x \\ J_y \end{pmatrix} = \begin{pmatrix} J_x \\ J_y \end{pmatrix} (J_x^2 + J_y^2). \qquad (26)$$

The other eigenvalues thus correspond to the square of the norm of the gradient of the image. Pixels where the gradient is zero (bad speckles) are also associated with zero eigenvalues.

M is also rank deficient, the dimension of its kernel is 2, a basis of its null space is well known:

$$\ker(M) = \operatorname{span} \begin{pmatrix} \mathbf{1} & 0 \\ 0 & \mathbf{1} \end{pmatrix}, \qquad (27)$$

where **1** is the vector filled with 1: the kernel of the scalar Laplace operator consists of constant functions. In fact a more efficient basis can be computed at a very low cost:

$$C = \begin{pmatrix} \frac{1}{\sqrt{s_{xx}}} & -\frac{s_{xy}s_b}{s_{xx}} \\ 0 & s_b\mathbf{1} \end{pmatrix} \text{ with } \begin{cases} s_{xx} = \mathbf{1}^T \mathbf{J}_x^2 \mathbf{1} \\ s_{xy} = \mathbf{1}^T \mathbf{J}_x \mathbf{J}_y \mathbf{1} \\ s_{yy} = \mathbf{1}^T \mathbf{J}_y^2 \mathbf{1} \\ s_b = 1/\sqrt{s_{yy} - s_{xy}^2/s_{xx}} \end{cases} \qquad (28)$$

It has the advantage to make the matrix $(C^T(A + \lambda M)C) = (C^T A C) = I$ for any $\lambda$.

## B Inverse of Laplacian on a rectangle with Neumann boundary condition

It is well known that plane waves $x \mapsto e^{i\omega \cdot x}$, with $\omega \in \mathbb{R}^2$, form a set of eigenfunctions for the Laplace operator in $\mathbb{R}^2$ with eigenvalues $-\|\omega\|^2$ (using the Euclidean norm). This can be equivalently formulated by saying that the Fourier transform diagonalizes the Laplacian. Hence, the powerful solution technique (in that case $\omega$ is the variable in the Fourier domain):

$$\begin{gathered} \Delta u + f = 0 \quad \text{in } \mathbb{R}^2 \\ u = \mathcal{F}^{-1} \left( \frac{\mathcal{F}(f)}{\|\omega\|^2} \right) \end{gathered} \qquad (29)$$

What is remarkable is that the eigenvectors are preserved by discretization. For instance, if we consider the classical 5-point stencil on a unit grid:

$$(\Delta_h u)(x_1, x_2) = u(x_1 + 1, x_2) + u(x_1 - 1, x_2) + u(x_1, x_2 + 1) + u(x_1, x_2 - 1) - 4u(x_1, x_2) \quad (30)$$

and one can check that

$$(\Delta_h e^{i\omega \cdot x})(x_1, x_2) = \underbrace{e^{i(x_1\omega_1 + x_2\omega_2)}}_{e^{i\omega \cdot x}} (e^{i\omega_1} + e^{-i\omega_1} + e^{i\omega_2} + e^{-i\omega_2} - 4) \qquad (31)$$

Now, considering a rectangular domain, the boundedness of the domain and the boundary conditions lead to selecting only certain eigenvalues, and eigenvectors are made out of a good combination of plane waves. Consider the unit square $[0, 1]^2$, the eigenvalues $\lambda_{n,m}$ and eigenvectors $v_{n,m}$ of the Laplacian with (homogeneous) Neumann boundary conditions are given by:

$$\begin{align} v_{n,m}(k,l) &= \cos\left(\frac{ml\pi}{M}\right) \cos\left(\frac{nk\pi}{N}\right) \tag{32} \\ \lambda_{n,m} &= 2\left(1 - \cos\left(\frac{n\pi}{N}\right)\right) + 2\left(1 - \cos\left(\frac{m\pi}{M}\right)\right) \nonumber \end{align}$$

As eigenvectors are cosine functions, the specialization of the Fourier transform to this case takes the name of discrete cosine transform (DCT).
The other eigenvalues thus correspond to the square of the norm of the gradient of the image. Pixels where the gradient is zero (bad speckles) are also associated with zero eigenvalues.

M is also rank deficient, the dimension of its kernel is 2, a basis of its null space is well known:

$$ \ker(\mathbf{M}) = \operatorname{span}\left(\begin{array}{cc} \mathbf{1} & 0 \\ 0 & \mathbf{1} \end{array}\right), \qquad (39) $$

where 1 is the vector filled with 1: the kernel of the scalar Laplace operator consists of constant functions. In fact a more efficient basis can be computed at a very low cost:

$$ \mathbf{C}_0 = \left( \begin{array}{cc} \frac{1}{\sqrt{s_{xx}}} \mathbf{1} & -\frac{s_{xy}s_b}{s_{xx}} \mathbf{1} \\ 0 & s_b \mathbf{1} \end{array} \right) \text{ with } \left\{ \begin{array}{l} s_{xx} = \mathbf{1}^T \mathbf{J}_x^2 \mathbf{1} = \sum_i j_{x,ii}^2 \\ s_{xy} = \mathbf{1}^T \mathbf{J}_x \mathbf{J}_y \mathbf{1} = \sum_i j_{x,ii} j_{y,ii} \\ s_{yy} = \mathbf{1}^T \mathbf{J}_y^2 \mathbf{1} = \sum_i j_{y,ii}^2 \\ s_b = 1/\sqrt{s_{yy} - s_{xy}^2/s_{xx}} \end{array} \right. \qquad (40) $$

This matrix is used as an augmentation in order to make sure that we only work in a space orthogonal to the kernel of M. It has the advantage to make the matrix $(\mathbf{C}_0^T(\mathbf{A} + \lambda\mathbf{M})\mathbf{C}_0) = (\mathbf{C}_0^T\mathbf{A}\mathbf{C}_0) = \mathbf{I}$ for any $\lambda$.

## 6.3 Quality of the preconditioner

We first wish to verify that preconditioning by regularization actually leads to better enforcement of the regularity. In Table 1, we can qualitatively compare the classical Jacobi approach of preconditioning by the diagonal of the operator $\mathbf{M}_{jac}^{-1} = \operatorname{diag}(\mathbf{A}_{\lambda})^{-1}$ and the proposed preconditioning by regularization. The increased regularity is particularly visible for low weight $\lambda$ and low precision $\varepsilon$ of the linear solver.

Preconditioning by the regularization operator thus makes it possible to make meaningful computations with low weight in the regularization and to solve with less precision, hence with fewer iterations. Nevertheless, one has to mention that our preconditioner is computationally more expensive per iteration than the diagonal one.

An interesting scenario unfolds. The regularization preconditioner promotes low frequency corrections. Indeed, it is associated with a fully populated matrix (never actually computed) and the search directions have naturally large wavelength. As iterations progress, higher frequency modes emerge introducing more and more details and irregularity. On the contrary, the Jacobi preconditioner is diagonal, and it naturally encourages (independent) details, only iterations make it possible to reveal the structure between neighboring pixels.

## 6.4 Ritz filtering

We analyze the solving process for high ($\lambda = 1000$) and low ($\lambda = 10$) levels of regularization. We use the second stopping criterion of Equation (12) with $\varepsilon = 10^{-5}$, which corresponds to a rather high level of convergence. The identified strain field are given in Figure 9.

We analyze the convergence in terms of compromise between the decrease of the error and the increase of the norm of the gradient of the solution which stems from the oscillations in the identified fields. Figure 7 presents two L-curves associated with high and low regularization. We use the natural CG-norms, please note that the position of the 0-abscissa is conventional because $\|\mathbf{x}_0 - \mathbf{x}\|_{\mathbf{A}}$ is unknown. The L-curves of the CG iterations (dotted lines) have similar shapes, like pieces of hyperbola. Due to the difference of magnitude, different scales had to be used: the error decreases four times less when the high regularization is used, and the norm of the solution remains 50 times smaller.

In order to better understand the convergence, we conduct a Ritz analysis. For $\lambda = 1000$, the convergence is attained in $m = 56$ iterations and as many Ritz vectors are computed. Table 2 presents a selection of these modes, sorted in decreasing order of Ritz value. The Ritz vectors resemble vibration modes with increasing number of anti-nodes. The first vectors are so regular that the hole is barely visible. The crack is only visible on the latest modes. This is bad (but logical) news because these are the most difficult modes to converge, thus they are probably bad approximations of actual eigenmodes, and the crucial mechanical information they carry is difficult to reuse.

We use formula (19) for the a posteriori filtering of the solution based on Ritz vectors. We present the L-curves in terms of modes included in the reconstruction. We show the curves in terms of full error and only taking into account the image error $(\mathbf{v}_i^T \mathbf{r}_{A,0})$ — they are almost overlaid on each other, a slight discrepancy only appears for high regularization. The shape of the Ritz L-curves corresponds to most of the modes (the highest) only slightly decreasing the error and almost not changing the norm, only the last modes, which contain
One just needs to take some care of the eigenvalue $\lambda_{0,0} = 0$, associated with the constant eigenvector. The classical solution is to work on functions with zero mean value and nullify the constant term in the transformed function.

We give the python code for the inverse of the discrete Laplacian on a rectangle with Neumann boundary conditions. This discrete Laplace operator can be directly invoked by the `laplace()` function from `scipy.ndimage` with default arguments (`border='reflect'`).

```python
import numpy as np
```

```python
from scipy.fftpack import dctn, idctn
```

```python
# Prepare transform of Laplacian for (N,M) images
mwx = 2 * (np.cos(np.pi*np.arange(0,N)/N)-1)
mwy = 2 * (np.cos(np.pi*np.arange(0,M)/M)-1)
[MWX, MWY] = np.meshgrid(mwx, mwy, indexing='ij')
MW = MWX + MWY
MW[0,0] = 1.
iMW = 1. / MW
iMW[0,0] = 0
```

```python
def SolveLaplaceNeumann(U,iMW): # U must have zero mean value
dctU = dctn(U, norm='ortho')
uhat = dctU * iMW
return(idctn(uhat,norm='ortho'))
```
<table><thead><tr><th>λ</th><th>ε</th><th>Diagonal Prec.</th><th>Regularization Prec.</th></tr></thead><tbody><tr><td>low</td><td>low</td><td></td><td></td></tr><tr><td>low</td><td>high</td><td></td><td></td></tr><tr><td>high</td><td>low</td><td></td><td></td></tr><tr><td>high</td><td>high</td><td></td><td></td></tr></tbody></table>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Table 1: $\epsilon_{xx}$ strain field (range = mean value ± 3 st.dev.). Comparison of the effect of preconditioning by diagonal (simple approach) vs by regularization, for different weights $\lambda \in \{1, 1000\}$ and linear solver precision $\epsilon \in \{10^{-2}, 10^{-3}\}$.</figcaption></figure>

the crack information, are associated with significant decrease of the error (but of course at the cost of much increased solution norm). More or less, if a corner was to be selected it would correspond to just suppressing the contribution of the last mode.

In order to better understand this behavior, we first analyze the convergence of the Ritz values by comparing the spectrum obtained at the last iteration with the one obtained just one iteration before, like was done in [14] in the case of a well posed problem. It appears that the largest Ritz values were quite well approximated and only the lowest part of the spectrum evolves (remember the Ritz values correspond to the inverse of the slope of the segments in the L-curve). In other words, even though the last iterations seem not to modify the solution much (accumulation of the dots in the upper left part on the CG L-curves), they play an important role in terms of estimation of the lower part of the spectrum, without adding lots of small eigenvalues.

To support this analysis, we conduct a Picard's study on Figure 8 which shows the distribution of the Ritz values ($\theta_i^{(m)}$) as well as the decomposition of the right-hand side on the eigenspace ($v_i^{(m)T} r_{A,0}$) and $\lambda(v_i^{(m)T} r_{M,0})$. It is worth recalling that low and highly regularized systems have the same spectrum, except that it is more sampled for the low regularization which requires two times more iterations to converge. The Ritz values are slowly decreasing and only the last 10% really decay, the low regularization is not associated with an overpopulation of the lowest part of the spectrum. What stands out is the fact that the right-hand side contributes almost equally on all modes (at least it does not decrease for larger Ritz values). Picard's theory
<figure><img src="image_15.png" /><figcaption>Table 2: $\epsilon_{xx}$ strain field for 8 Ritz vectors out of the 56 computed ($\epsilon = 10^{-5}$, $\lambda = 1000$).</figcaption></figure>

thus suggests that we should stop the reconstruction when the Ritz values start to decay. This is not possible in our case since the crack is mostly represented in this part of the spectrum.

By the way, Figure 8 permits to compare the smallest Ritz value $\theta_{\min}$ with the regularization parameter $\lambda$. The case that we called “low regularization” corresponds to $\lambda$ being negligible with respect to the small Ritz value $\theta_m^{(m)}$, and thus only marginally modifying the active Ritz spectrum. On the contrary, the high regularization corresponds to a $\lambda > \theta_m^{(m)}$ which means that the lower part of the spectrum of $A_\lambda$ is flattened relative to that of $A$.

<figure><img src="image_16.png" /><figcaption>Figure 7: L-curves (same $\epsilon = 10^{-5}$) compared with Ritz post-treatment, for different regularization intensity $\lambda$. Note that different scales are used on the plots.</figcaption></figure>

## 6.5 Subspace recycling

Even though it appears that Ritz filtering is difficult to apply to the studied system, we can still benefit from Ritz vectors to accelerate the solution. As a sequence of linear systems with identical matrix has to be solved, it is natural to augment the system with the previously generated Ritz vectors by concatenating $C \leftarrow (C V_m)$. Indeed, augmentation comes with optimized block operations that make augmenting by one vector much cheaper than one iteration in particular using (30) and (31).

<table><caption>Table 3: Performance of recycling</caption><tr><td>Aug.</td><td>0</td><td>10</td><td>20</td><td>30</td><td>40</td><td>50</td><td>60</td><td>70</td><td>max (77)</td></tr><tr><td>Iter.</td><td>77</td><td>64</td><td>57</td><td>49</td><td>44</td><td>41</td><td>40</td><td>40</td><td>38</td></tr><tr><td>Time (s)</td><td>11.7</td><td>10.1</td><td>8.3</td><td>7.1</td><td>6.6</td><td>6.3</td><td>6.2</td><td>6.3</td><td>6.6</td></tr></table>
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>Figure 8: Spectral analysis of the system for $\epsilon = 10^{-5}$ and different regularization intensity $\lambda$. A 5-width median filter was used to smooth out the contribution curves.</figcaption></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>Figure 9: Identified $\epsilon_{xx}$ field, for different regularization intensity $\lambda$, with $\epsilon = 10^{-5}$.</figcaption></figure>

Table 3 illustrates the performance of recycling for the nonlinear system to be solved on the full image at the end of identification. The first linear system, only augmented by the kernel of the preconditioner is solved in 77 iterations. Then a certain portion of the Ritz vectors is used to augment the next 8 linear systems (same matrix, different right-hand sides). As the augmentation results in an excellent initialization, we use a criterion in terms of absolute value of $\|\mathbf{r}_j\|_{\mathbf{M}}$ to halt the iterations because other comparison as given in Equation (12) might use an unfair reference. We measure the performance in terms of gain in iterations, and in computational time (measures are conducted on a upper mid-range laptop with Nvidia RTXA2000 graphic card). The gain in terms of iterations is moderate, with best obtained for small augmentation space (at most 1.3 iterations per augmentation vector, for 10 vectors). In terms of time, the optimal is obtained for augmentation space of 80%-90% of available vectors, with a global CPU time divided by almost 2 (this time includes all the extra cost associated with computing and using Ritz vectors). This size of subspace agrees with what we observed on the stability of the largest Ritz values in the L-curves plots.

## 6.6 Tuning of $\lambda$

It is often hard to automatize the selection of the regularization intensity $\lambda$. Picard's plots like in Figure 8 permit to put $\lambda$ in relation with the spectrum of the preconditioned operator and thus to understand the effect of the regularization in terms of flattened spectrum. Still, the final judge is often the expert's impression of a strain map, and it is convenient to compute maps associated with several ($\lambda_p$) at low cost.

As presented in (32), after solving one linearized system regularized by $\lambda_0$, we post-process the solution for any $\lambda_p$ at the simple cost of computing the associated right-hand side $b_{\lambda_p}^k(\mathbf{x}_{\lambda_p}^k)$ (which depends on the history of the nonlinear solution for $\lambda_p$), and basic linear algebra operations.

Figure 10 presents the solution deduced for $\lambda_1 = 1$ from initial computations with different $\lambda_0$ (in $\{10, 1000, 10000\}$)
<figure><img src="image_21.png" /><figcaption>Figure 10: Costless postprocessing for different ($\lambda_i$): top, initial computation with $\lambda_0 \in \{10, 100, 1000\}$; bottom, direct solution with $\lambda = 1$ and solutions deduced with $\lambda_1 = 1$. $\epsilon_{xx}$ strain field.</figcaption></figure>

and $\varepsilon = 10^{-4}$. Again, a median filter was applied after the Ritz reconstruction. It seems that the Ritz vectors make it possible to postprocess a reasonable solution with up to $\lambda_1 \simeq \lambda_0/1000$. The reconstructed strain field appears to be much less smooth than the original computation (with $\lambda_0$) while less noisy than the direct low-regularization computation with $\lambda = 1$. If the deduced solution is not fully satisfying, it can still be used as an excellent initialization for a regular computation.

# 7 Conclusion

In this paper, we have studied how preconditioning and Tikhonov regularization could be efficiently combined in an augmented preconditioned conjugate gradient. We have shown that this association makes sense from a physical point of view and it made it possible to combine filtering, recycling of subspaces, and postprocessing of all regularized solutions at negligible cost. This gives a favorable framework to apply criteria like the L-curve or Picard's analysis.

First, the solver was applied to the boundary completion problem. It was verified that regularization was a good preconditioner and that Ritz analysis made it possible to approximate pretty efficiently the L-curve and to build extremely informative Picard plots.

The solver was then applied to a problem of the optical flow reconstruction which introduced the extra difficulty of nonlinearity and the fact that the most important information was buried in the lower part of the spectrum. Preconditioning by the regularization made it possible to start from a regular estimation and enhance details through iterations. Satisfying results were obtained on actual measurements from digital image correlation of a mechanical test. Postprocessed solutions at negligible cost were still relevant for regularization weight $\lambda$ divided by up to 1000.

An obvious next step for this work is to consider inexact preconditioners, that is to say when the $M^{-1}$ matrix in the preconditioning step of the algorithm is only an approximation of the inverse of the regularization matrix $M$ in the operator. This would make the method applicable on a much broader class of problems.

# References

[1] Owe Axelsson and Igor Kaporin. Error norm estimation and stopping criteria in preconditioned conjugate gradient iterations. *Numerical Linear Algebra with Applications*, 8(4):265-286, 2001.

[2] Owe Axelsson and Gunhild Lindskog. On the rate of convergence of the preconditioned conjugate gradient method. *Numerische Mathematik*, 48:499-523, 1986.

[3] Faker Ben Belgacem. Why is the Cauchy problem severely ill-posed? *Inverse Problems*, 23(2):823-836, 2007.

[4] Faker Ben Belgacem and Henda El Fekih. On Cauchy's problem: I. A variational Steklov-Poincaré theory. *Inverse Problems*, 21(6):1915, 2005.
[5] Gilles Besnard, François Hild, and Stéphane Roux. “Finite-Element” displacement fields analysis from digital images: application to Portevin–Le Châtelier bands. *Experimental Mechanics*, 46(6):789–803, 2006.

[6] Gilles Besnard, François Hild, and Stéphane Roux. Finite-element displacement fields analysis from digital images: application to portevin-le châtelier bands. *Experimental Mechanics*, 46:789–804, 2006.

[7] Alessandro Buccini, Marco Donatelli, and Lothar Reichel. Iterated Tikhonov regularization with a general penalty term. *Numerical Linear Algebra with Applications*, 24(4):e2089, 2017. e2089 nla.2089.

[8] Ahmed Chabib, Jean-François Witz, Pierre Gosselet, and Vincent Magnier. GCPU OpticalFlow: a GPU accelerated python software for strain measurement. *SoftwareX*, 26:101688, March 2023.

[9] Ahmed Chabib, Jean-François Witz, Pierre Gosselet, and Vincent Magnier. The impact of metrics in mechanical imaging. *Strain*, 61:e1249, 2025.

[10] Andrew Chapman and Youssef Saad. Deflated and augmented Krylov subspace techniques. *Numerical Linear Algebra with Applications*, 4(1):43–66, 1997.

[11] David Chin-Lung Fong and Michael Saunders. CG versus MINRES: An empirical comparison. *Sultan Qaboos University Journal for Science*, 17(1):44–62, 2012.

[12] Zdeněk Dostál. Conjugate gradient method with preconditioning by projector. *International Journal of Computer Mathematics*, 23:315–323, 1988.

[13] Renaud Ferrier, Mohamed L. Kadri, and Pierre Gosselet. The Steklov-Poincaré technique for data completion: Preconditioning and filtering. *International Journal for Numerical Methods in Engineering*, 116(4):270–286, 2018.

[14] Pierre Gosselet, Christian Rey, and Julien Pebrel. Total and selective reuse of Krylov subspaces for the resolution of sequences of nonlinear structural problems. *International Journal for Numerical Methods in Engineering*, 94(1):60–83, 2013.

[15] Per Christian Hansen. The truncated SVD as a method for regularization. *BIT Numerical Mathematics*, 27(4):534–553, 1987.

[16] Per Christian Hansen. The discrete Picard condition for discrete ill-posed problems. *BIT Numerical Mathematics*, 30(4):658–672, 1990.

[17] Per Christian Hansen. Analysis of discrete ill-posed problems by means of the L-curve. *SIAM review*, 34(4):561–580, 1992.

[18] Magnus R. Hestenes and Eduard Stiefel. Methods of conjugate gradients for solving linear systems. *Journal of research of the national bureau of standards*, 49(6):409–436, 1952.

[19] Zhongxiao Jia and G.W. Stewart. On the convergence of the Ritz values, Ritz vectors and refined Ritz vectors. Technical Report 3896, Institute of Advanced Computer Studies, Department of Computer Science, University of Maryland at College Park, 1999.

[20] Mohamed Larbi Kadri, Jalel Ben Abdallah, and Thouraya Nouri Baranger. Identification of internal cracks in a three-dimensional solid body via Steklov-Poincaré approaches. *Comptes Rendus Mécanique*, 339(10):674–681, 2011.

[21] Louis Kovalevsky and Pierre Gosselet. A quasi-optimal coarse problem and an augmented Krylov solver for the Variational Theory of Complex Rays. *International Journal for Numerical Methods in Engineering*, 2015.

[22] Vladimir Arkad'evich Kozlov, Vladimir Gilelevich Maz'ya, and AV Fomin. An iterative method for solving the Cauchy problem for elliptic equations. *Zhurnal Vychislitel'noi Matematiki i Matematicheskoi Fiziki*, 31(1):64–74, 1991.

[23] Gérard Meurant, Jan. Papež, and Petr P. Tichý. Accurate error estimation in CG. *Numerical Algorithms*, 88:1337–1359, 2021.

[24] Vladimir Alekseevich Morozov. The error principle in the solution of operational equations by the regularization method. *Zhurnal Vychislitel'noi Matematiki i Matematicheskoi Fiziki*, 8(2):295–309, 1968.

[25] Andreas Neubauer. Ill-posed problems and the conjugate gradient method: Optimal convergence rates in the presence of discretization and modelling errors. *Journal of Inverse and Ill-posed Problems*, 30(6):905–915, 2022.
[26] Stanley Osher, Martin Burger, Donald Goldfarb, Jinjun Xu, and Wotao Yin. An iterative regularization method for total variation-based image restoration. *Multiscale Modeling & Simulation*, 4(2):460–489, 2005.

[27] Jean-Charles Passieux and Robin Bouclier. Classic and inverse compositional Gauss-Newton in global DIC. *International Journal for Numerical Methods in Engineering*, 119(6):453–468, 2019.

[28] Yousef Saad. *Iterative methods for sparse linear systems*. SIAM, 2003.

[29] Yousef Saad. *Numerical Methods for Large Eigenvalue Problems*, volume 66 of *Classics in Applied Mathematics*. SIAM, Philadelphia, USA, revised edition, 2011.

[30] Deqing Sun, Stefan Roth, and Michael J. Black. Secrets of optical flow estimation and their principles. In *2010 IEEE Computer Society Conference on Computer Vision and Pattern Recognition*, pages 2432–2439, 2010.

[31] Andreři Nikolaevich Tikhonov and Vasilić Iakovlevich Arsenin. *Solutions of ill-posed problems*. Wiley, 1977.

[32] Wenli Wang, Gangrong Qu, Caiqin Song, Youran Ge, and Yuhan Liu. Tikhonov regularization with conjugate gradient least squares method for large-scale discrete ill-posed problem in image restoration. *Applied Numerical Mathematics*, 204:147–161, 2024.

[33] Xiangtuan Xiong, Xuemin Xue, and Zhi Qian. A modified iterative regularization method for ill-posed problems. *Applied Numerical Mathematics*, 122:108–128, 2017.

# A Inverse of Laplacian on a rectangle with Neumann boundary condition

It is well known that plane waves $x \mapsto e^{i\omega \cdot x}$, with $\omega \in \mathbb{R}^2$, form a set of eigenfunctions for the Laplace operator in $\mathbb{R}^2$ with eigenvalues $-\|\omega\|^2$ (using the Euclidean norm). This can be equivalently formulated by saying that the Fourier transform diagonalizes the Laplacian. Hence, the powerful solution technique (in that case $\omega$ is the variable in the Fourier domain):

$$ 
\begin{array}{l}
\Delta u + f = 0 \quad \text{in } \mathbb{R}^2 \\
u = \mathcal{F}^{-1} \left( \frac{\mathcal{F}(f)}{\|\omega\|^2} \right)
\end{array}
\qquad (41) $$

What is remarkable is that the eigenvectors are preserved by discretization. For instance, if we consider the classical 5-point stencil on a unit grid:

$$ (\Delta_h u)(x_1, x_2) = u(x_1 + 1, x_2) + u(x_1 - 1, x_2) + u(x_1, x_2 + 1) + u(x_1, x_2 - 1) - 4u(x_1, x_2) \quad (42) $$

and one can check that

$$ (\Delta_h e^{i\omega \cdot x})(x_1, x_2) = \underbrace{e^{i(x_1\omega_1 + x_2\omega_2)}}_{e^{i\omega \cdot x}} (e^{i\omega_1} + e^{-i\omega_1} + e^{i\omega_2} + e^{-i\omega_2} - 4). \quad (43) $$

Now, considering a rectangular domain, the boundedness of the domain and the boundary conditions lead to selecting only certain eigenvalues, and eigenvectors are made out of a good combination of plane waves. Consider the unit square $[0, 1]^2$, the eigenvalues $\lambda_{n,m}$ and eigenvectors $v_{n,m}$ of the Laplacian with (homogeneous) Neumann boundary conditions are given by:

$$ 
\begin{aligned}
v_{n,m}(k,l) &= \cos\left(\frac{ml\pi}{M}\right) \cos\left(\frac{nk\pi}{N}\right) \\
\lambda_{n,m} &= 2\left(1 - \cos\left(\frac{n\pi}{N}\right)\right) + 2\left(1 - \cos\left(\frac{m\pi}{M}\right)\right)
\end{aligned}
\qquad (44) $$

As eigenvectors are cosine functions, the specialization of the Fourier transform to this case takes the name of discrete cosine transform (DCT).

One just needs to take some care of the eigenvalue $\lambda_{0,0} = 0$, associated with the constant eigenvector. The classical solution is to work on functions with zero mean value and nullify the constant term in the transformed function.

Figure 11 gives the python code for the inverse of the discrete Laplacian on a rectangle with Neumann boundary conditions. This discrete Laplace operator can be directly invoked by the `laplace()` function from `scipy.ndimage` with default arguments (`border='reflect'`).
```python
import numpy as np
from scipyURRE import dctn, idctn

# Prepare transform of Laplacian for (N,M) images
mwx = 2 * (np.cos(np.pi*np.arange(0,N)/N)-1)
mwy = 2 * (np.cos(np.pi*np.arange(0,M)/M)-1)
[MWX, MWY] = np.meshgrid(mwx, mwy, indexing='ij')
MW = MWX + MWY
MW[0,0] = 1.
iMW = 1. / MW
iMW[0,0] = 0

def SolveLaplaceNeumann(U, iMW):
    # U must have zero mean value
dctU = dctn(U, norm='ortho')
uhat = dctU * iMW
return(idctn(uhat, norm='ortho'))
```

<caption>Figure 11: Inverse discrete Laplacian on a rectangle (Neumann bcs) in python</caption>