Interferometric phase linking (IPL) has become a prominent technique for processing images of areas containing distributed scaterrers in SAR interferometry.

Traditionally, IPL consists in estimating consistent phase differences between all pairs of SAR images in a time series from the sample covariance matrix of pixel patches on a sliding window.

This paper reformulates this task as a covariance fitting problem: in this setup, IPL appears as a form of projection of an input covariance matrix so that it satisfies the phase closure property.

Given this modular formulation, we propose an overview of covariance matrix estimates, regularization options, and matrix distances, that can be of interest when processing multi-temporal SAR data.

In particular, we will observe that most of the existing IPL algorithms appear as special instances of this framework.

We then present tools to efficiently solve related optimization problems on the torus of phase-only complex vectors: majorization-minimization and Riemannian optimization.

We conclude by illustrating the merits of different options on a real-world case study.

§ INTRODUCTION

Recent remote sensing missions (Sentinel-1, UAVSAR, TerraSAR-X, etc.) have brought an unprecedented amount of available synthetic aperture radar (SAR) images time series.

For interferometric SAR (InSAR), these systematic and regular acquisitions enabled the utilization of multi-temporal techniques (MT-InSAR), which significantly enhanced the accuracy of Earth displacement estimation.

Under the assumption of distributed targets, interferometric phase-linking (IPL) has emerged as a fundamental methodology to estimate phase differences from all available SAR acquisitions within the dataset [1, 2, 3, 4, 5, 6, 7, 8, 9, 10].

The driving idea behind this technique is to leverage the redundancy of the time series
in order to compensate for the coherence loss between images over time, as illustrated in Fig. <ref>.

The initial formulation of IPL was obtained from the perspective of approximate maximum likelihood estimation assuming a complex circular Gaussian model.

This estimation procedure then appears as an optimization problem aiming to retrieve phase estimates from the sample covariance matrix of each pixel patch.

Subsequent works motivated many variants, e.g., by improving coherence pre-estimation step [5], using compression schemes [7], relaxing the optimization problem [6], or assuming non-Gaussian models [11, 12].

An overview of advances in this scope is presented in [13].

In essence, all IPL algorithms aim to recover a property called phase closure (referring to the continuity of phase differences over the time series) within the sample covariance matrix.

In this paper, we leverage this reinterpretation to reformulate IPL as a covariance fitting problem [14, 15, 16, 17, 18, 19, 20, 21]: the task is thus expressed as fitting constrained phases to the modulus of any plug-in estimate of the covariance matrix.

This approach, referred to as COFI-PL, offers a concise and modular formulation.

It encompasses the maximum likelihood approaches as a special case, as well as many other generalizations concerning the construction of the covariance matrix plug-in estimate, and the fitting objective function.

In this scope, we overview staple building blocks from the state-of-the-art and relate them to existing IPL algorithms.

Furthermore, we consider new options bringing promising results on a real-world case study, notably the use of least-squares fitting objective, and the phase-only sample correlation matrix.

As all of the considered methods lead to the construction of an optimization problem over the torus of phase-only complex vectors $\mathbb{T}_p$, we further investigate two optimization frameworks suited to this space.

Contrarily to the majority of existing works that focus on optimizing the phases directly, this offers several advantages: it simplifies tedious computations of trigonometric functions to simpler matrix operations, and enables us to fully harness the geometric structure of the constrained space (for example, invariance of the phases modulo $2\pi$ is not an appearing issue when considering $\mathbb{T}_p$ directly).

The first framework is majorization-minimization [22, 23] on $\mathbb{T}_p$.

Its use is limited to quadratic problems on $\mathbb{T}_p$, which interestingly includes classical maximum likelihood [24, 25] and least-squares based [9] IPL formulations.

The main appeal of this framework lies in its ability to yield cost-efficient algorithms (here similar to modified power-methods), with guaranteed monotonicity, and with no required hyperparameter such as a step-size selection.

Hence, they appear suited to large-scale implementation.

The second framework is Riemannian optimization on $\mathbb{T}_p$ [26, 27, 28].

It is more flexible as it allows to work with any objective function.

It also opens the path to acceleration schemes and generalizations from many algorithms of the Riemannian optimization framework [29, 30].

The rest of the paper is organized as follows:

Section <ref> presents IPL and the corresponding covariance matrix structure.

Section <ref> presents the generic construction of a covariance fitting IPL (COFI-PL) optimization problem.

This framework involves three construction modules: some design options suited to SAR data are respectively presented in Sections <ref> (covariance matrix estimation),
Section <ref> (covariance matrix regularization), and Section <ref> (matrix distance objective function).

Section <ref> and Section <ref> respectively introduce the majorization-minimization and Riemannian optimization on $ \mathbb{T}_p$ to solve
the resulting COFI-PL problems.

Section <ref> overviews existing IPL algorithms and links them to COFI-PL formulations.

Section <ref> presents an application of the proposed methods to a real-world case study.

Naive interferogram formed between two SAR images (left) and interferogram formed with interferometric phase linking by leveraging an interleaved time series between those two images (right).

§ BACKGROUND

§.§ Notations

The following convention is adopted: italic indicates a scalar quantity, lower case boldface indicates a vector quantity and upper case boldface a matrix.

For a complex-valued number $x = a+ib=re^{i\theta}\in\mathbb{C}$,
$\mathfrak{Re}(x)=a$, $\mathfrak{Im}(x)=b$, ${\rm mod}(x)=r$, ${\rm arg}(x)=\theta$ and $\phi_{\mathbb{T}}(x)= e^{i\theta}$.

These operators overload to matrix entries by being applied element-wise.

$\mathcal{H}_p^{++}$ (resp. $\mathcal{S}_p^{++}$) is the set of hermitian (resp. symmetric) positive definite matrices.

$\mathbb{T}_p$ is the torus of phase only complex vectors of dimension $p$ as in (<ref>).

The operators ${\rm tr}(\cdot)$ and $|\cdot|$ return respectively the trace and the determinant of a matrix.

The entrywise complex conjugation is denoted $\cdot^*$, while the transpose (resp. transpose conjugate) operation is denoted $\cdot^\top$ (resp. $\cdot^H$).

The Hadamard product is denoted $\circ$.

The matrix $\mathbf{I}$ denotes the identity matrix of appropriate dimension.

A circular multivariate Gaussian vector of mean $\boldsymbol{\mu}$ and covariance matrix $\mathbf{\Sigma}$ is denoted $\mathbf{x}\sim \mathcal{CN}(\boldsymbol{\mu},\mathbf{\Sigma})$.

§.§ InSAR covariance matrix structure

Stack of $p$ co-registered SAR images and multi-looking window: gray pixels represent the current local patch, denoted $\{\mathbf{x}_i\}_{i=1}^n$.

From a given stack of $p$ co-registered SAR images, we consider a sliding window that processes all local patches of $n$ multivariate pixels, as illustrated in Figure <ref>.

A pixel patch is denoted as $\{\mathbf{x}_i\}_{i=1}^n$, with $\mathbf{x}_i \in \mathbb{C}^{p},~\forall i\in [\![1,n]\!]$.

A single multivariate pixel $\mathbf{x}_i$ contains a time-series (in chronological order) of the $p$ snapshots, i.e.

\begin{equation}
\mathbf{x}_i = \left[~ x_i^1 ,~\cdots,~x_i^{p}  ~ \right]^{\top} ~\in \mathbb{C}^p.
\end{equation}

The set $\{\mathbf{x}_i\}_{i=1}^n$ is assumed to be a homogeneous patch containing $n$ adjacent pixels with similar scattering and statistical properties.

From the standard physical considerations about distributed scatterers in SAR, we have that the first and second order moments as follows:

\begin{equation}\label{eq:moments}
\begin{aligned}
& \mathbb{E} \left[ x^q \right] = 0,~\forall~k\in [\![1,p]\!] \\
& \mathbb{E} \left[ x^{q}  ( x^{\ell} )^* \right] =
\upsilon_{q,\ell}
\sigma_q \sigma_{\ell} e^{j(\theta_q - \theta_{\ell})}
,~\forall (q,\ell) \in [\![1,p]\!]^2 \\
\end{aligned}
\end{equation}

$\bullet$ $\sigma_q^2 = \mathbb{E}\left[ x^q(x^q)^H\right] \in \mathbb{R}^+$ is the variance of $x^q$.
The corresponding vector of standard deviations is denoted $\boldsymbol{\sigma} = \left[ \sigma_1 ,~ \cdots,~ \sigma_{p}  \right]$.

$\bullet$ $\upsilon_{q,\ell}  \in [0 , 1]$
is the coherence coefficient between $x^q$ and $x^{\ell}$.
The corresponding coherence matrix is denoted $\boldsymbol{\Upsilon}$, with entries $[ \boldsymbol{\Upsilon} ]_{q,\ell} = \upsilon_{q,\ell} \in [0,1]$. Also remark that $[ \boldsymbol{\Upsilon} ]_{\ell,\ell}=1, \forall l\in [\![1,p]\!]$.

$\bullet$ $\theta_q$ is the phase instant $q$.
We denote the phase vector $\boldsymbol{\theta} = \left[ \theta_{1},~\cdots,~\theta_{p}\right]$, and the corresponding vector of complex phases[
Because phase-only complex vectors [26] will be extensively used in this work, we need to distinguish two objects related to the polar decomposition of complex numbers.

For $x = a+ib=re^{i\theta}\in\mathbb{C}$, $\theta$ is indifferently called the phase, angle, or argument.
However, $e^{i\theta}$ will be specifically referred to as the complex phase.] is

\begin{equation}\label{eq:w_vector}
\mathbf{w}_{\boldsymbol{\theta}} =
\left[ e^{j\theta_{1}},~\cdots,~e^{j\theta_{p}}\right] ~\in \mathbb{T}_p,
\end{equation}

\begin{equation} \label{eq:phaseonlytorus}
\mathbb{T}_p = \{ \mathbf{w}\in\mathbb{C}^p ~|~|[\mathbf{w}]_q| =1, ~\forall q\in [\![1,p]\!] \}
\end{equation}

is the torus of phase-only complex vectors.

By convention, we will use the reference $\theta_1 = 0$, which is equivalent to $[\mathbf{w}_{\boldsymbol{\theta}}]_1=1$.

The covariance structure in (<ref>) is expressed in matrix form as

\begin{equation} \label{eq:cov_struct_1}
\mathbb{E} \left[ \mathbf{x}  \mathbf{x}^H \right]
= \mathbf{\Sigma}
= {\rm diag}(\mathbf{w}_{\boldsymbol{\theta}}) {{\boldsymbol{\Psi}}}
{\rm diag}(\mathbf{w}_{\boldsymbol{\theta}})^H
= \boldsymbol{\Psi}  \circ (\mathbf{w}_{\boldsymbol{\theta}}
\mathbf{w}_{\boldsymbol{\theta}}^H),
\end{equation}

where $\boldsymbol{\Psi} =    {\boldsymbol{\Upsilon}} \circ\boldsymbol{\sigma} \boldsymbol{\sigma}^\top$ is the coherence matrix scaled by the variance coefficients.
We can also notice that this decomposition coincides with the modulus-argument decomposition, i.e.:

\begin{equation}\label{eq:mod_arg_decomp}
\mathbf{\Sigma} = {\rm mod}(\mathbf{\Sigma})  \circ \phi_\mathbb{T} (\mathbf{\Sigma}) \overset{\Delta}{=} \boldsymbol{\Psi}  \circ (\mathbf{w}_{\boldsymbol{\theta}}
\mathbf{w}_{\boldsymbol{\theta}}^H).
\end{equation}

where we used the complex phase extraction operator defined as $\phi_\mathbb{T} : x {=} r e^{i\theta} \mapsto e^{i\theta}$ (which extends naturally to matrices by being applied entry-wise).

From the covariance matrix expression in (<ref>) and (<ref>), we observe that $\mathbf{\Sigma}$
cannot be any covariance matrix in $\mathcal{H}_p^{++}$, as it exhibits the particular phase structure $ {\rm arg}(\mathbf{\Sigma}) = \mathbf{w}_{\boldsymbol{\theta}} \mathbf{w}_{\boldsymbol{\theta}}^H$.

Indeed, if we denote the phase differences between two images indexed $q$ and $\ell$ as $\Delta_{q,\ell} = \theta_q - \theta_{\ell}$, we have

\begin{equation}
\Delta_{q,\ell} +
\Delta_{\ell,j} +
\Delta_{j,q} = 0
\end{equation}

that is satisfied for all triplet $\{q,\ell,j\}$.

Such relationship translates directly into $\mathbf{\Sigma}$ as

\begin{equation} \label{eq:phase_closure}
{\rm arg}  (\mathbf{\Sigma}_{q\ell})   +
{\rm arg}  (\mathbf{\Sigma}_{\ell j})  +
{\rm arg}  (\mathbf{\Sigma}_{j q}) = 0
\end{equation}

because ${\rm arg} (\mathbf{\Sigma}_{q\ell}) = \Delta_{q,\ell}$ from (<ref>).

The aforementioned property is referred to as phase closure, or phase consistency.
It is an important property in MT-InSAR, as it is related to the continuity of physical phenomena, such as Earth displacement.

§.§ Interferometric phase linking

In practice, the true covariance matrix $\mathbf{\Sigma}$ of the data is unknown.

The interferometric phases have thus to be estimated solely from the sample set $\{\mathbf{x}_i\}_{i=1}^n$.

From equation (<ref>), a naive approach would consists in computing the sample covariance matrix

\begin{equation} \label{eq:SCM}
\mathbf{S} = \frac{1}{n} \sum_{i=1}^n \mathbf{x}_i \mathbf{x}_i^H ,
\end{equation}

and identifying the phase difference from the arguments of its entries, i.e.,

\begin{equation} \label{eq:naive_ipl}
\hat{\Delta}_{q,\ell} =     {\rm arg}  \ (\mathbf{S}_{q\ell})
\end{equation}

Unfortunately, this simple estimate is relatively inaccurate (especially when $|q-\ell|$ increases, due to the temporal decorrelation).

It also provides a series of phase difference estimates that do not satisfy the phase closure (<ref>).

More accurate estimation procedures consist rather in directly estimating the vector of complex phases $\mathbf{w}_{\boldsymbol{\theta}}$ (or equivalently, the phase vector $\boldsymbol{\theta}$) from $\mathbf{S}$ by leveraging the prior structure (<ref>).

The process is referred to as phase triangulation, or interferometric phase linking (IPL), for which numerous algorithms have been developed over the years (cf. [13] for a recent overview).

In this scope, the next section proposes a general framework capable of encompassing most of the existing methods, providing a modular structure for their extensions.

§ COFI-PL FRAMEWORK

This section reformulates IPL as a generic covariance fitting problem, whose corresponding framework will be referred to as COFI-PL.

Covariance fitting (or covariance matching) is a widely employed technique in array processing, which consists in refining the structure of an input covariance matrix estimator by minimizing a projection criterion [14, 15, 16, 17, 18, 19, 20, 21].

In this setup, IPL can be interpreted as projecting (according to some distance or divergence) the input estimate to the set of matrices that satisfy the phase closure property.

In particular, we will see that most of the established MLE-inspired algorithms (e.g., from [1, 2, 4]) appear as a special case of COFI-PL when considering the Kullback-Leibler (KL) divergence as a fitting criterion.

Given any (possibly regularized) plug-in estimate of the covariance matrix $\mathbf{\Sigma}$, denoted $\tilde{\mathbf{\Sigma}}$, that does not satisfy the phase closure, our objective is to determine the “best” phase-constrained fitting given some distance criterion linking $\tilde{\mathbf{\Sigma}}$ and its modulus denoted

\begin{equation}
\tilde{\mathbf{\Psi}} \overset{\Delta}{=} {\rm mod}(\tilde{\mathbf{\Sigma}}).
\end{equation}

The problem is formulated as

\begin{equation}\label{eq:COFI-PL}
\begin{array}{c l}
\underset{\mathbf{w}_{\boldsymbol{\theta}}}{\rm minimize}
d^2 (  \tilde{\mathbf{\Sigma}}
\tilde{\boldsymbol{\Psi}}
\circ
\mathbf{w}_{\boldsymbol{\theta}}\mathbf{w}_{\boldsymbol{\theta}}^H
\\
{\rm subject~to}
\mathbf{w}_{\boldsymbol{\theta}} \in \mathbb{T}_p
\\
\theta_1=0
\end{array}
\end{equation}

where $d$ is a matrix distance (or divergence) that will be specified later on.

Note that the only variable is the vector $\mathbf{w}_{\boldsymbol{\theta}}$, so we will use the compact notation

\begin{equation}
f^d_{\tilde{\mathbf{\Sigma}}} ( \mathbf{w}_{\boldsymbol{\theta}}) =
d^2 (  \tilde{\mathbf{\Sigma}}
\tilde{\boldsymbol{\Psi}}
\circ
\mathbf{w}_{\boldsymbol{\theta}}\mathbf{w}_{\boldsymbol{\theta}}^H
\end{equation}

for the objective in (<ref>).

Also remark that the objective function in (<ref>) is invariant to a constant phase-shift of all entries in $\mathbf{w}_{\boldsymbol{\theta}}$, thus the constraint $\theta_1=0$ can be discarded, and achieved a posteriori by subtracting $\theta_1$ to all the optimized phases.

The generic COFI-PL problem is finally expressed as

\begin{equation}\label{eq:COFI-PL2}
\begin{array}{c l}
\underset{\mathbf{w}_{\boldsymbol{\theta}}}{\rm minimize}
f^d_{\tilde{\mathbf{\Sigma}}} ( \mathbf{w}_{\boldsymbol{\theta}})
\\
{\rm subject~to}
\mathbf{w}_{\boldsymbol{\theta}} \in \mathbb{T}_p.
\end{array}
\end{equation}

This formulation then offers a multitude of options concerning:

* The construction of the covariance matrix plug-in $\hat{\mathbf{\Sigma}}$
and its possible regularization $\tilde{\mathbf{\Sigma}}$.

* The choice of the matrix distance $d$ according to a geometry of interest.

* The optimization method to address the constraint $\mathbf{w}_{\boldsymbol{\theta}} \in \mathbb{T}_p$ efficiently.

With these options specified, the corresponding instance of COFI-PL then defines a complete processing chain from the data $\{\mathbf{x}_i\}_{i=1}^n$ to the phase estimates, as illustrated in Figure <ref>.

The following sections <ref>, <ref>, and <ref> will provide an overview of practical design options suited to SAR data.

Then sections <ref> and <ref> will present two optimization frameworks that can be used to solve (<ref>) under the constraint $\mathbf{w}_{\boldsymbol{\theta}} \in \mathbb{T}_p$.

A modular COFI-PL chain of process: the covariance matrix plug-in $\tilde{\mathbf{\Sigma}}$
is estimated from the samples $\{\mathbf{x}_i\}_{i=1}^n$ with a possible regularization.
The choice of a matrix distance $d$ then specifies the covariance fitting optimization problem (<ref>).
Interferometric phase estimates are obtained by solving this problem using optimization methods on the $n$-torus of phase-only vectors $\mathbb{T}_p$.

§ UNSTRUCTURED COVARIANCE MATRIX ESTIMATION

The first module requires to construct an estimate of the covariance matrix $\hat{\mathbf{\Sigma}}$ from the sample set $\{\mathbf{x}_i\}_{i=1}^n$.

In this setup, it is interesting to link the chosen estimation method to underlying assumptions on the statistical model of the data [31].

Notably, SAR data can be non-Gaussian [32, 25], which motivates the use of robust estimation methods.

This section overviews relevant plug-in estimates for SAR data, and their corresponding model assumptions.

§.§ Sample covariance matrix

A common assumption in SAR is to consider that each sample is independent and identically distributed according to a centered complex circular Gaussian model, denoted $\mathbf{x}\sim\mathcal{CN}(\mathbf{0},\mathbf{\Sigma})$.

This model is particularly relevant for low resolution SAR images, i.e., when we can assume that each pixel gathers the sum of the contributions from many scatterers, and apply the central limit theorem.

In this case, the pixel patch has the negative log-likelihood

\begin{equation} \label{eq:gauss_likelihood}
\mathcal{L}_{\mathcal{N}}(\{\mathbf{x}_i\}_{i=1}^n|\mathbf{\Sigma})  \propto {\rm tr} (\mathbf{S}\mathbf{\Sigma}^{-1} ) + \log | \mathbf{\Sigma} |
\end{equation}

with $\mathbf{S}$ as in (<ref>).

In this case, the maximum likelihood estimator of $\mathbf{\Sigma}$ is the sample covariance matrix $\mathbf{S}$.

This estimate is the most widely used plug-in in IPL [2, 4].

§.§ $M$-estimators of the scatter

At high resolution, the Gaussian assumption is often no longer valid, and empirical histograms tend to exhibit heavy-tails (see, e.g., [32, 25]).

The centered circular complex elliptically symmetric (CES) distributions provides a better fit for such data.

CES models, denoted $\mathbf{x}\sim\mathcal{CES}(\mathbf{0},\mathbf{\Sigma},g)$, correspond to the following negative log-likelihood function:

\begin{equation} \label{eq:ces_likelihood}
\mathcal{L}_{\mathcal{E}}^g
\propto
\frac{1}{n}\sum_{i=1} \rho (   \mathbf{x}_i^H \mathbf{\Sigma}^{-1} \mathbf{x}_i    )
+ \log | \mathbf{\Sigma} |
\end{equation}

with $\rho(t) = - \log g(t)$, and where $g$ is referred to as the density generator.

The special case $g(t)=e^{-t}$ yields back the Gaussian model, while many other options allow for modelling other heavy-tailed distributions [33].

An $M$-estimator of the scatter $ \mathbf{\Sigma}_{M}$ is then defined as the solution of the following fixed-point equation

\begin{equation}
\label{eq:mest}
\mathbf{\Sigma}_{M} = \frac{1}{n}\sum_{i=1}^{n} u(\mathbf{x}_i^H  \mathbf{\Sigma}_{M}^{-1}\mathbf{x}_i)\mathbf{x}_i\mathbf{x}_i^H
\end{equation}

where $u$ is a real-valued weight function on $[0,\infty)$.

Such solution exists and is unique under conditions on $u$ and the sample support (notably, $n>p$) discussed in [34, 33].

When $u(t) = - g'(t) / g (t)$, (<ref>) is the maximum likelihood estimator of $\mathbf{\Sigma}$ for $\mathbf{x} \sim \mathcal{CES}(\mathbf{0},\mathbf{\Sigma}, g)$.

Otherwise, $M$-estimators are robust to a mismatch, and we can use a function $u$ that is not necessarily linked to $g$.

Most notably, Tyler's estimator [35], based on the function $u_{T}(t) = p/t$, is distribution free, in the sense that its distribution does not depend on the underlying density generator $g$.

In practice, $M$-estimators still offer better estimation performance compared to the sample covariance matrix as long as the function $u$ assesses for potentially heavy-tailed distribution.

Specifically for IPL, these estimators have been considered as plug-ins in [11, 12].

§.§ Correlation-based plug-in

First, recall that ${\rm mod} (\mathbf{\Sigma})=\boldsymbol{\Psi} =  {\boldsymbol{\Upsilon}} \circ\boldsymbol{\sigma} \boldsymbol{\sigma}^\top$, where ${\boldsymbol{\Upsilon}}$ is the coherence matrix and $\boldsymbol{\sigma}$ is the vector of standard deviations.

Remark that the scaling by $\boldsymbol{\sigma}$ does not impact the phase structure of $\mathbf{\Sigma}$ in (<ref>).

This can motivate the use of a correlation (rather than covariance) matrix estimate $\hat{\mathbf{C}}$ as plug-in, in order to mitigate issues related to the amplitude fluctuation in SAR images (unbalanced backscattered power among all the images).

Such matrix is built from any covariance matrix estimate $\hat{\mathbf{\Sigma}}$ as

\begin{equation} \label{eq:C_standardization}
\hat{\mathbf{C}} = {\rm diag}(\hat{\mathbf{\Sigma}})^{-1/2} \hat{\mathbf{\Sigma}}\: {\rm diag}(\hat{\mathbf{\Sigma}})^{-1/2}.
\end{equation}

Most notably, the sample correlation matrix, i.e., (<ref>) built with $\hat{\mathbf{\Sigma}}=\mathbf{S}$, was successfully leveraged in [4], and other works that implicitly standardize the data before computing $\mathbf{S}$.

Still, this estimator can be sensitive to non-Gaussian distributions.

A solution is to turn to robust estimators of the correlation matrix, for which many options exist [36, 37].

We will focus on a simple and highly robust one, that is obtained from the phase-only sample correlation matrix

\begin{equation}
\mathbf{T} = \frac{1}{n} \sum_{i=1}^n \mathbf{y}_i \mathbf{y}_i^H,
\label{eq:phase_only_SCM}
\end{equation}

with $\mathbf{y}=\phi (\mathbf{x})$.

As this estimator projects each entry of the samples to the unit sphere, it is inherently robust to any underlying distributions of the modulus of the marginals $x^q,~\forall q\in[\![1,p]\!]$ (cf. [38] for an example of such distributions).

§ COVARIANCE MATRIX REGULARIZATION

The second module consists in applying a regularization process to the estimator provided by the first one.

Indeed, the sample support $n$ is limited by the size of the sliding window (that sets the spatial resolution).

For long time series, this can lead to situations where $n\simeq p$ or $n<p$, in which the plain estimate $\hat{\mathbf{\Sigma}}$ is inaccurate.

Applying some form of regularization to the estimate in this case greatly improves the accuracy of phase estimation.

This section thus overviews how to construct a regularized estimator $\tilde{\mathbf{\Sigma}}$ from a plug-in estimate $\hat{\mathbf{\Sigma}}$ with different approaches motivated by IPL, and assumptions on SAR data.

These come at the cost of regularization parameter selection, so we also discuss references that address this issue.

§.§ Shrinkage to identity

At low sample support, covariance matrix estimates are usually ill conditioned ($n\simeq p$), or even not invertible ($n<p$).

This poses a major issue in IPL because most of the fitting distances are constructed from the inverse of the plug-in estimate (cf. Section  <ref>).

A practical solution is to operate a shrinkage of the estimate to a scaled identity matrix, in order to leverage some form of bias-variance trade-off.

A popular formulation of such regularization preserves the scale (i.e., trace) of the estimate, and is defined for any plug-in estimate $\hat{\mathbf{\Sigma}}$ as

\begin{equation} \label{eq:LWshrink}
\tilde{\mathbf{\Sigma}}(\beta)
\beta \hat{\mathbf{\Sigma}}
(1-\beta) \frac{{\rm tr}(\hat{\mathbf{\Sigma}}) }{p} \mathbf{I}
\end{equation}

with $\beta\in [0,1]$.

This regularization is also referred to as spectral shrinkage, as it shrinks eigenvalues of the plug-in estimate towards their mean.

When the plug-in estimate is the sample covariance matrix $\mathbf{S}$, the adaptive selection of the regularization parameter $\beta$ for minimizing the mean squared error has been studied in [39] (assuming finite $4^{\text{th}}$ order moments), [40] (assuming Gaussian data), and [41] (assuming elliptically distributed data).

Extension to the $M$-estimators as in (<ref>) has been studied in [42].

§.§ Low-rank approximation

The empirical spectrum of multivariate SAR pixel patches often exhibits a low-rank structure (see e.g. experiments in [43, 25]).

This means that most of the data variance lies in a rank-$k$ linear subspace, which justifies the use of principal component analysis[
In practice the rank $k$ is usually set fixed for the whole image, i.e. for processing all pixel patches indifferently.
Though it comes at a heavy computational cost, this process could be refined by using adaptive rank estimation, e.g., using model order selection methods [44].
] [45].

The corresponding decomposition of the covariance matrix has often been
used in IPL to improve the quality of the plug-in estimate [4].

Let the eigenvalue decomposition of this plug-in estimate be denoted $\hat{\mathbf{\Sigma}} \overset{\rm EVD}{=} \sum_{r=1}^p \hat{\lambda}_r \hat{\mathbf{u}}_r \hat{\mathbf{u}}_r^H$.

The shrinkage approaches that leverage a low-rank structure within this decomposition consider either the rank-$k$ approximation of the plug-in estimate

\begin{equation} \label{eq:low_rank_approx}
\tilde{\mathcal{P}}_k (\hat{\mathbf{\Sigma}} ) = \sum_{r=1}^k \hat{\lambda}_r \hat{\mathbf{u}}_r \hat{\mathbf{u}}_r^H ,
\end{equation}

or its projection on the set of rank-$k$ plus scaled identity, defined as

\begin{equation}  \label{eq:low_rank_approx2}
{\mathcal{P}}_k (\hat{\mathbf{\Sigma}} ) =
\sum_{r=1}^k \hat{\lambda}_r \hat{\mathbf{u}}_r \hat{\mathbf{u}}_r^H
\sum_{r=k+1}^p \bar{\lambda}_r \hat{\mathbf{u}}_r \hat{\mathbf{u}}_r^H ,
\end{equation}

with $\bar{\lambda}_r = \sum_{r=k+1}^p \hat{\lambda}_r /(p-k)$.

This second option is advocated in this work for two main reasons.

The first is theoretical, as applying $\mathcal{P}_k$ to the sample covariance matrix $\mathbf{S}$ yields the maximum likelihood estimator of a low-rank Gaussian signal plus white Gaussian noise [46, 47].

The second is practical: the operator $\mathcal{P}_k$ gives an invertible matrix and improves the conditioning of the plug-in estimate by uplifting its lowest eigenvalue to the average of the $p-k$ lowest ones.

Similar to the shrinkage to identity, this property is instrumental to compute distances $d$ whose expressions involve matrix inverses.

§.§ Covariance matrix tapering

SAR image stacks suffer from targets decorrelation over time, which is why IPL exploits the redundancy brought by all image pairs in order to construct the interferograms.

However, time frames that exhibit a prohibitively low coherence should intuitively be excluded from this construction, which motivated the development of methods to determine which pairs of images are exploited or disregarded [48, 49].

In this scope, a simple idea consists in processing only pairs contained in a sliding temporal window of bandwidth $b$ (where the coherence is assumed to remain high enough).

Within the IPL framework, this translates into forcing a banded structure in the covariance matrix plug-in.

Let $\mathbf{W}(b)$ be a banding-type[Note that other tapering templates
could be envisioned [50, 51], but we focus only on the most relevant to the considered applications.] tapering matrix with bandwidth $b$ [52, 53]

\begin{equation}
[\mathbf{W}(b)]_{ij} =
\left\{
\begin{array}{cl}
1  & \text{if}~|i-j|\leq b\\
0  & \text{otherwise}.
\end{array}
\right.
\end{equation}

Covariance matrix tapering, also referred to as Hadamard regularization, involves producing the regularized estimate

\begin{equation} \label{eq:tapering}
\tilde{\mathbf{\Sigma}} (b) = \mathbf{W}(b) \circ \hat{\mathbf{\Sigma}},
\end{equation}

which then naturally exhibits a banded structure.

Adaptive procedures for the bandwidth selection can be found in [52, 53, 54].

The joint use of shrinkage to identity and tapering was advocated for IPL in [8].

The optimal (in the sense of mean-squared error) adaptive selection of the parameters $(\beta,b)$ for this regularization was studied in [54] and applied to IPL in [10].

§.§ Joint estimation and regularization

For completeness, we mention that some covariance matrix plug-ins can be constructed by “merging” the estimation and regularization modules.

A main example is that $M$-estimators do not exist for $n<p$, so the regularization cannot be performed afterwards, and has to be included within the robust estimation process.

Regularized $M$-estimators, i.e., minimizers of $\mathcal{L}_{\mathcal{E}}^g
(\{\mathbf{x}_i\}_{i=1}^n|\mathbf{\Sigma})$ in (<ref>)
plus an additive penalty term, were studied in [55, 56, 57, 58], and their use for IPL was discussed in [59].

Low-rank structured $M$-estimators, i.e., minimizers of $\mathcal{L}_{\mathcal{E}}^g
(\{\mathbf{x}_i\}_{i=1}^n|\mathbf{\Sigma})$ in (<ref>) under the low-rank structure constraint as in [46, 47], were investigated [60, 43].

§ MATRIX DISTANCES

The previous sections presented two modules dedicated to the construction of a (possibly regularized) plug-in estimate of the covariance matrix $\tilde{\mathbf{\Sigma}}$.

From this estimate, the third module formulates the optimization problem (<ref>), which established a covariance matching type IPL.

This process simply requires to select a matrix distance, for which many options are available.

This section overviews prominent ones, and explores their connections to statistical assumptions and existing IPL methods.

§.§ Kullback-Leibler divergence and maximum likelihood estimation approaches

The Kullback-Leibler (KL) divergence measures the dissimilarity between two probability density functions.

Its expression between two centered Gaussian distributions $\mathbf{x} \sim \mathcal{CN}(0, \mathbf{\Sigma}_1)$ and $\mathbf{x}\sim \mathcal{CN}(0, \mathbf{\Sigma}_2)$ is:

\begin{multline}
\label{eq:KL_gauss}
{\rm KL}(\mathcal{CN}(0, \mathbf{\Sigma}_1 )~||~\mathcal{CN}(0, \mathbf{\Sigma}_2  ))
\\ = {\rm tr}(\mathbf{\Sigma}_2^{-1}\mathbf{\Sigma}_1) + \log |\mathbf{\Sigma}_2\mathbf{\Sigma}_1^{-1}|-p,
\end{multline}

which provides a matrix divergence between $\mathbf{\Sigma}_1$ and $\mathbf{\Sigma}_2$.

Setting $\mathbf{\Sigma}_1 = \tilde{\mathbf{\Sigma}} $, $\mathbf{\Sigma}_2 = \tilde{\boldsymbol{\Psi}}
\circ \mathbf{w}_{\boldsymbol{\theta}}  \mathbf{w}_{\boldsymbol{\theta}}^H $, keeping $\mathbf{w}_{\boldsymbol{\theta}}$ as the only variable, and simplifying the expression, we obtain the KL-IPL objective function as

\begin{equation} \label{eq:kl_COFI-PL}
f^{\rm KL}_{\tilde{\mathbf{\Sigma}}}( \mathbf{w}_{\boldsymbol{\theta}})  = \mathbf{w}_{\boldsymbol{\theta}}^H
( \tilde{\boldsymbol{\Psi}}^{-1} \circ \tilde{\mathbf{\Sigma}} )
\mathbf{w}_{\boldsymbol{\theta}}.
\end{equation}

This objective function is the most widely employed for IPL, especially when choosing the sample covariance matrix as plug-in estimate.

However, the expression (<ref>) with $\tilde{\mathbf{\Sigma}} = \mathbf{S}$ is usually obtained from the perspective of maximum likelihood estimation under the Gaussian model when assuming that ${\rm mod}(\mathbf{\Sigma})$ (coherence matrix scaled by the variance coefficients) is known [1, 2, 4].

Owing to similar expressions of the Gaussian log-likelihood in (<ref>) and the KL divergence (<ref>) when using $\mathbf{S}$ as the plug-in estimate, the two approaches fall back on the same objective in this case.

There is a slight difference between these approaches when using regularization, which is discussed further in the overview Section <ref>.

A last remark is that the KL divergence is not symmetric, so it is not a proper matrix distance.
Its symmetric counterpart (i.e., inverting the roles of $\mathbf{\Sigma}_1$ and $\mathbf{\Sigma}_2$), and its symmetrized version could also be envisioned as options.

Still, we focus only on the chosen formulation because of its direct link to the Gaussian maximum likelihood estimator, and well-established IPL algorithms.

§.§ Frobenius norm and least-squares estimator

The euclidean distance, also referred to as “flat metric”, between two symmetric matrices is defined as

\begin{equation}
d_{E}^2 ( \mathbf{\Sigma}_1, \mathbf{\Sigma}_2) = || \mathbf{\Sigma}_1 - \mathbf{\Sigma}_2 ||^2_F .
\end{equation}

By setting $\mathbf{\Sigma}_1 = \tilde{\mathbf{\Sigma}} $, $\mathbf{\Sigma}_2 = \tilde{\boldsymbol{\Psi}}
\circ \mathbf{w}_{\boldsymbol{\theta}}  \mathbf{w}_{\boldsymbol{\theta}}^H $, while keeping $\mathbf{w}_{\boldsymbol{\theta}}$ as the only variable, and simplifying the expression, we obtain the LS-IPL (for least squares) objective function as

\begin{equation} \label{eq:ls_COFI-PL}
f^{\rm LS}_{\tilde{\mathbf{\Sigma}}} (\mathbf{w}_{\boldsymbol{\theta}} )
-2 \mathbf{w}_{\boldsymbol{\theta}}^H
( \tilde{\boldsymbol{\Psi}} \circ \tilde{\mathbf{\Sigma}} )
\mathbf{w}_{\boldsymbol{\theta}}
+ {\rm const.},
\end{equation}

This objective function appears as similar to the KL-IPL problem (<ref>), i.e., a quadratic form on $\mathbb{T}_p$.

However, KL-IPL is formulated as a minimization of the quadratic form involving $ \tilde{\boldsymbol{\Psi}}^{-1} \circ \tilde{\mathbf{\Sigma}} $, while LS-IPL is formulated as a maximization (because of the minus sign in (<ref>)) one involving $ \tilde{\boldsymbol{\Psi}} \circ \tilde{\mathbf{\Sigma}} $.

Though the LS-IPL formulation has driven more research due to its maximum likelihood grounds [13], recent studies [9, 10] evidenced the practical use of $f^{\rm LS}$.

A notable interest of this objective function is that no covariance matrix inversion is required: besides the reduction of the computational complexity, it can also mitigate inaccuracies arising from the poor conditioning of the plug-in estimate (without requiring any regularization process).

§.§ Weighted Frobenius norm and EXIP criterion

The weighted Frobenius norm between two symmetric matrices is defined as:

\begin{equation}
d_{\rm WF}^2
\mathbf{\mathbf{H}}^{-1/2}
\mathbf{\mathbf{H}}^{-1/2}
\end{equation}

where $\mathbf{H}$ is a whitening-type weight matrix to be fixed.

Among possible choices, we will focus on $\mathbf{H}=\tilde{\mathbf{\Sigma}}$,
that yields the extended invariance principle (EXIP) and COMET-type estimators, which hold interesting statistical properties [14, 15, 16, 17, 18, 19, 20, 21].

Such EXIP approach has been used in InSAR in [61, 62, 63] (though not always directly related to IPL formulations).

For this choice, setting $\mathbf{\Sigma}_1 = \tilde{\mathbf{\Sigma}} $ and $\mathbf{\Sigma}_2 = \tilde{\boldsymbol{\Psi}}
\circ \mathbf{w}_{\boldsymbol{\theta}}  \mathbf{w}_{\boldsymbol{\theta}}^H $ gives the following WLS-IPL (for weighted least squares) objective function

\begin{equation}
\begin{aligned}
f^{\rm WLS}_{\tilde{\mathbf{\Sigma}}}
& =
\mathbf{I}
\tilde{\mathbf{\Sigma}}^{-1/2}
\circ \mathbf{w}_{\boldsymbol{\theta}} \mathbf{w}_{\boldsymbol{\theta}}^H)
\tilde{\mathbf{\Sigma}}^{-1/2}
\end{aligned}
\end{equation}

which unfortunately, cannot be simplified into a simpler quadratic form as KL-IPL and LS-IPL.

§ MAJORIZATION-MINIMIZATION ON $\MATHBB{T}_P$

The options discussed in the previous section allow us to construct optimization problems as in (<ref>) in order to perform IPL.

This section presents majorization-minimization on $\mathbb{T}_p$ in order to address the resolution of these problems.

This optimization framework has been successfully leveraged for KL-IPL [25, 24] and LS-IPL [9].

Though it is not applicable to all objective functions, its major practical interest lies in the fact that it leads to simple and scalable algorithms.

§.§ Majorization-minimization

Majorization-minimization is briefly reviewed here with notations that match the problem (<ref>) for convenience.

More details on this framework can be found in [22, 23].

We consider an optimization problem of the form

\begin{equation}
\begin{array}{c l}\label{eq:problem_phases2}
\underset{{\mathbf{w}}\in\mathbb{T}_p}{\rm minimize}
f(\mathbf{w}) .
\end{array}
\end{equation}

The majorization-minimization algorithm is an iterative optimization procedure that operates with two steps:

* Majorization: at current point $\mathbf{w}_t$, find a surrogate function $g(\cdot |\mathbf{w}_t)$ so that it is tangent to the objective, $f(\mathbf{w}_t) = g(\mathbf{w}_t | \mathbf{w}_t)$, and majorizes it, i.e.,

\begin{equation}
f(\mathbf{w})\leq  g(\mathbf{w} | \mathbf{w}_t),~ \forall \mathbf{w} \in \mathbb{T}_p
\end{equation}

* Minimization: obtain the next iterate as

\begin{equation}
\mathbf{w}_{t+1} = {\rm argmin}_{\mathbf{w}\in \mathbb{T}_p } ~ g(\mathbf{w} | \mathbf{w}_t).
\end{equation}

This algorithm enjoys nice convergence properties [64] (being constrained to the compact set $\mathbb{T}_p$ can be accounted for with the same arguments as in [65, 66]).

It notably ensures a monotonic decrement of the objective function at each step.

The main interest of this approach is that it can yield a sequence of sub-problems that are easily solved if the surrogate functions are well constructed.

§.§ Lemmas for quadratic forms of phase only complex vectors

This section presents useful surrogates functions and their closed form minimizer for quadratic forms on $\mathbb{T}_p$.

Let $\mathbf{H}\succcurlyeq \mathbf{0}$ be a hermitian positive semi-definite matrix, we have the following lemmas.

The concave quadratic form $f:\mathbf{w} \mapsto - \mathbf{w}^H \mathbf{H} \mathbf{w}$ is majorized at point $\mathbf{w}_t$ by the surrogate function

\begin{equation}
g ( \mathbf{w} | \mathbf{w}_t )  =
- 2 \mathfrak{Re}
\{ \mathbf{w}^H
\underbrace{\mathbf{H} \mathbf{w}_t }_{\tilde{\mathbf{w}}_t}
\}
+ {\rm const. }
\end{equation}

with equality at point $\mathbf{w}_t$.

A concave function lies below its tangent curves, so it can be majorized at any point by its first order Taylor expansion.

The convex quadratic form $f:\mathbf{w} \mapsto \mathbf{w}^H \mathbf{H} \mathbf{w}$ is majorized on $\mathbb{T}_p$ at point $\mathbf{w}_t$ by the surrogate function

\begin{equation}
g ( \mathbf{w} | \mathbf{w}_t )  =
2 \mathfrak{Re}
\{ \mathbf{w}^H
\underbrace{(\mathbf{H} - \lambda_{\rm max}^{\mathbf{H}} \mathbf{I} ) \mathbf{w}_t }_{-\tilde{\mathbf{w}}_t}
\}
+ {\rm const. }
\label{eq:linmajofquadcvx}
\end{equation}

where $\lambda_{\rm max}$ be the largest eigenvalue of $\mathbf{H}$, and with equality achieved at $\mathbf{w}_t$.

We first notice that, when $\mathbf{w}\in\mathbb{T}_p$, we have the relation

\begin{equation} \label{eq:cvx_ccv_relation}
\mathbf{w}^H (\mathbf{H} - \lambda_{\rm max}^{\mathbf{H}} \mathbf{I} ) \mathbf{w}
\mathbf{w}^H \mathbf{H} \mathbf{w} -
\underbrace{p \lambda_{\rm max}^{\mathbf{H}}}_{ {\rm const.}}.
\end{equation}

Hence, the objective function restricted to $\mathbb{T}_p$ coincides with the expression of a concave quadratic form up to a constant.
Optimizing either side of the equality in (<ref>) over $\mathbb{T}_p$ thus yields the same solution.
The quadratic form $\mathbf{w}^H (\mathbf{H} - \lambda_{\rm max}^{\mathbf{H}} \mathbf{I} ) \mathbf{w} $ is concave, thus it can be majorized at point $\mathbf{w}_t$ by its first order Taylor expansion as in (<ref>).

The two previous lemmas show that any quadratic form can be majorized on $\mathbb{T}_p$ by linear functions.

For these linear surrogates functions, the minimization step on $\mathbb{T}_p$ can be solved in closed form thanks to the following lemma.

The solution to the minimization problem

\begin{equation}
\begin{array}{c l}\label{eq:linpr_torus}
\underset{{\mathbf{w}}\in\mathbb{T}_p}{\rm minimize}
- \mathfrak{Re} \left\{ \mathbf{w}^H \bar{\mathbf{w}}_t \right\}
\end{array}
\end{equation}

is obtained as

\begin{equation}
\mathbf{w}^\star = \phi_\mathbb{T}( \bar{\mathbf{w}}_t )
\end{equation}

with the operator $\phi_\mathbb{T} : x {=} r e^{i\theta} \mapsto e^{i\theta}$ (that
extends naturally to matrices by being applied entry-wise).

The problem requires maximizing the sum of $p$ independent inner products in $\mathbb{C}$ of the form $\langle {e}^{i\theta_j}, [\bar{\mathbf{w}}_t]_j \rangle$, which is solved in closed-from by aligning the phase of each entry as $\theta_j = {\rm arg} ([\tilde{\mathbf{w}}_t]_j) $.

§.§ Application to KL-IPL and LS-IPL

We can remark that LS-IPL and KL-IPL are formulated as minimization problems over $\mathbb{T}_p$ whose objective functions are quadratic forms, i.e., as instances of a generic problem

\begin{equation} \label{eq:generif_torus_qf}
\begin{array}{c l}
\underset{{\mathbf{w}}\in\mathbb{T}_p}{\rm minimize}
\mathbf{w}^H \mathbf{M} \mathbf{w}
\end{array}
\end{equation}

* KL-IPL involves the matrix $\mathbf{M}_{\rm KL}= \tilde{\mathbf{\Psi}}^{-1} \circ \tilde{\mathbf{\Sigma}} $, that is generally positive semi-definite (cf. remark below).

Thus we can use Lemma <ref> and Lemma <ref> in order to obtain the algorithm summarized in the box Algorithm <ref>.

* LS-IPL involves the matrix $\mathbf{M}_{\rm LS}= - \tilde{\mathbf{\Psi}} \circ \tilde{\mathbf{\Sigma}}$, that is negative semi-definite (cf. remark below).

Thus we can use Lemma <ref> and Lemma <ref> in order to obtain the algorithm summarized in the box Algorithm <ref>.

The Hadamard product of two positive semi-definite matrices is semi-definite.

However, the modulus of a positive semi-definite matrix is not guaranteed to stay positive semi-definite (some counterexamples can be found numerically).

It means that $\mathbf{M}_{\rm KL}$ (resp. $\mathbf{M}_{\rm LS}$) is not necessarily positive semi-definite (resp. negative semi-definite) by construction.

Though these corner cases were not actually observed in practice, it can be interesting to implement some safeguards checks.

As remedy, it is always possible to shift the eigenvalues to ensure positiveness, as done in Lemma <ref>.

The other option is to split the negative and positive part of the matrix, i.e., decomposing $\mathbf{M} = \mathbf{M}_+  + \mathbf{M}_-$, with $\mathbf{M}_+\succcurlyeq \mathbf{0}$ and $\mathbf{M}_- \preccurlyeq  \mathbf{0}$ by eigenvalue decomposition, before applying Lemma 1 and 2 separately to each resulting quadratic form.

Majorization-minimization for KL-IPL

[1]
Entry: $\tilde{\mathbf{\Sigma}} \in \mathbb{C}^{p\times p}$ (plug-in), $\mathbf{w}_1 \in \mathbb{T}_p$ (starting point)
Compute $\mathbf{M} = {\rm mod} (\tilde{\mathbf{\Sigma}} )^{-1} \circ \tilde{\mathbf{\Sigma}}  $  and $\lambda_{\rm max}^{\mathbf{M}}$
Compute ${\bar{\mathbf{w}}_t} = (\lambda_{\rm max}^{\mathbf{M}} \mathbf{I} - \mathbf{M} ) \mathbf{w}_t $
Update $\mathbf{w}_t = \phi_{\mathbb{T}} \{ \bar{\mathbf{w}}_t \} $
Output: $\hat{\mathbf{w}}_{\boldsymbol{\theta}} = \mathbf{w}_{\rm end} \in \mathbb{T}_p$

Majorization-minimization for LS-IPL

[1]
Entry: $\tilde{\mathbf{\Sigma}} \in \mathbb{C}^{p\times p}$ (plug-in), $\mathbf{w}_1 \in \mathbb{T}_p$ (starting point)
Compute $\mathbf{M} = {\rm mod}(\tilde{\mathbf{\Sigma}} ) \circ \tilde{\mathbf{\Sigma}} $
Compute ${\bar{\mathbf{w}}_t} =  \mathbf{M} \mathbf{w}_t $
Update $\mathbf{w}_t = \phi_{\mathbb{T}} \{ \bar{\mathbf{w}}_t \} $
Output: $\hat{\mathbf{w}}_{\boldsymbol{\theta}} = \mathbf{w}_{\rm end} \in \mathbb{T}_p$

§.§ EMI-type EVD relaxations

The structure of the problem in (<ref>) is reminiscent of the computation of an eigenvector.

Motivated by this analogy, the EMI algorithm [6] relaxed the constraint $\mathbf{w} \in \mathbb{T}_p$ into the following problem:

\begin{equation}   \label{eq:EMIrelax}
\begin{array}{c l}
\underset{{\mathbf{w}}}{\rm minimize}
\mathbf{w}^H \mathbf{M} \mathbf{w} \\
{\rm subject~to} & \mathbf{w}^H \mathbf{w} = 1 .
\end{array}
\end{equation}

The solution to this problem is the eigenvector associated with the lowest (resp. largest) eigenvalue of $\mathbf{M}$ if it is hermitian positive (resp., negative) definite.

Thus it can be easily computed, and the corresponding phase estimates are then identified directly from the complex phases of the entries of this eigenvector.

In general, this relaxed estimator is not as accurate as the one brought by the proper resolution of (<ref>).

However, it has been evidenced to be satisfactory in some practical cases [6, 10], justifying its use in terms of performance versus computational load trade-off.

§ RIEMANNIAN OPTIMIZATION ON $\MATHBB{T}_P$

Though majorization-minimization provides simple and elegant algorithms, its application is limited to IPL objective functions that are based on quadratic forms.

Because $\mathbb{T}_p$ (<ref>) is a smooth manifold, more general objective functions can be handled with the theory of Riemannian optimization [29, 30].

This framework has been leveraged for optimization on $\mathbb{T}_p$ in beamforming application [26, 27, 28]. However, to the best of our knowledge, it has not been considered within the context of IPL.

Its main interest compared to direct optimization of the argument vector $\boldsymbol{\theta}$ is that it greatly simplifies the calculations, while inherently accounting for the geometry and invariances of the considered space.

This section presents basic tools that allow developing first-order based methods such as the Riemannian counterparts of steepest descent, conjugate gradient, or Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithms.

The torus of phase-only complex vector $\mathbb{T}_p$ in (<ref>)
is a compact smooth manifold embedded in $\mathbb{C}^p$.

Its tangent space at point $\mathbf{w}$ is denoted $T_{\mathbf{w}} \mathbb{T}_p$, and defined

\begin{equation}
\label{eq:tangent_torus}
T_{\mathbf{w}} \mathbb{T}_p =
\left\{
\boldsymbol{\xi} \in \mathbb{C}^p
\mathfrak{Re}\{ \boldsymbol{\xi} \circ \mathbf{w}^* \} = \mathbf{0}
\right\}
\end{equation}

which can intuitively be recovered by identifying it to a product of tangent space of $p$ complex circles (representing each entry of the vector $\mathbf{w})$.

Then $\mathbb{T}_p$ can be turned into a Riemannian submanifold <cit.> by endowing each tangent space with the Euclidean metric

\begin{equation}
\label{eq:metric_torus}
\begin{array}{l c c l}
\langle \cdot , \cdot \rangle_\mathbf{w} : & T_{\mathbf{w}} \mathbb{T}_p \times T_{\mathbf{w}} \mathbb{T}_p &\rightarrow& \mathbb{R}  \\
& \boldsymbol{\xi}, \boldsymbol{\eta} & \mapsto & \mathfrak{Re}\{ \boldsymbol{\xi}^H \boldsymbol{\eta}  \}
\end{array}
\end{equation}

We then consider optimization problem of the form (<ref>), that has no obvious closed-form solution on $\mathbb{T}_p$.

In order to evaluate this solution, we resort to iterative
methods, i.e., methods that yield a sequence of iterates $\mathbf{w}_t\in \mathbb{T}_p$ from a starting point $\mathbf{w}_0 \in \mathbb{T}_p$.

This sequence is constructed so that it eventually converges to
a critical point of the objective in (24).

When the variable is constrained to lie in the manifold $\mathbb{T}_p$, first-order based Riemannian optimization methods operate as follows:

* At iterate $\mathbf{w}_t$, a descent direction in the tangent space, denoted $\boldsymbol{\xi}_t\in T_{\mathbf{w}} \mathbb{T}_p$, is computed by leveraging the Riemannian gradient of the objective function.

* The direction descent  $\boldsymbol{\xi}_t$ is used to obtain the next iterate $\mathbf{w}_{t+1}$. This is achieved through a retraction on $\mathbb{T}_p$, which is an operator that maps tangent vectors back onto the manifold.

This generic procedure is illustrated in Figure <ref>.

For the first step, the steepest direction descent is given by the Riemannian gradient, which is defined according to the metric (<ref>).

The Riemannian gradient of a function $f: \mathbb{T}_p \rightarrow \mathbb{R} $ at point $\mathbf{w}_t$ is the unique vector defined as

\begin{equation}
\langle {\rm grad} f(\mathbf{w}_t) , \boldsymbol{\xi} \rangle_{\mathbf{w}_t} =
{\rm D} f(\mathbf{w}_t ) \left[ \boldsymbol{\xi} \right],
\end{equation}

where ${\rm D} f(\mathbf{w}_t ) \left[ \boldsymbol{\xi} \right]$ is the directional derivative of $f$ with respect to $\mathbf{w}_t$ in direction $\boldsymbol{\xi} $.

This directly yields the Riemannian gradient as the orthogonal projection of the Euclidean gradient onto the tangent space, i.e.:

\begin{equation} \label{eq:riemgrad}
{\rm grad} f(\mathbf{w}_t)
\nabla f(\mathbf{w}_t) - \mathfrak{Re} \{ \nabla f(\mathbf{w}_t)^* \circ \mathbf{w}_t \} \circ \mathbf{w}_t,
\end{equation}

where $\nabla f(\mathbf{w}_t)$ is the Euclidean gradient of the objective
function at $\mathbf{w}_t$.

Momentum-based, or acceleration methods generally combine the gradient of several iterates.

Because these objects belong to different tangent spaces in our context, we need
to meaningfully map tangent vectors from one point to another.

This is achieved thanks to parallel transport <cit.>, which is obtained in our context with

\begin{equation}
\mathcal{T}_{\mathbf{w}_{t}\leftarrow\mathbf{w}_{t-1}}^{\mathbb{T}_p}
= \boldsymbol{\xi} - \mathfrak{Re}( \boldsymbol{\xi} \circ \mathbf{w}_t^* ) \circ \mathbf{w}_t.
\end{equation}

where $\mathcal{T}_{\mathbf{w}_{t}\leftarrow\mathbf{w}_{t-1}}^{\mathbb{T}_p}$ is used to denote the transport of a tangent vector $\boldsymbol{\xi}\in T_{\mathbf{w}_{t-1}} \mathbb{T}_p$ to the tangent space $T_{\mathbf{w}_{t}} \mathbb{T}_p$.

Given a direction descent $\boldsymbol{\xi}_t \in T\mathbf{w}_t\mathbb{T}_p $,
the second step can be performed by defining a retraction, i.e., an operator $\mathcal{R}_{\mathbf{w}_t}: T\mathbf{w}_t \mathbb{T}_p \to \mathbb{T}_p $ that maps tangent vectors back to the manifold.

In our case, the euclidean projection on $\mathbb{T}_p$ is a practical candidate, so we can use

\begin{equation} \label{eq:retraction}
\mathcal{R}_{\mathbf{w}_t} ({\boldsymbol{\xi}}_t) = \phi_\mathbb{T}(\mathbf{w}_t + \boldsymbol{\xi}_t ),
\end{equation}

where $\phi_\mathbb{T}$ was defined after (<ref>).

The gradient, parallel transport, and retraction are enough to derive most first-order based algorithms by transposing from standard optimization methods.

As an example, the Riemannian conjugate gradient is given in the box Algorithm <ref>, while the Euclidean gradient for covariance fitting cost computed with LS and KL can be obtained from the relations:

\begin{equation}
\nabla  \mathbf{w}^H \mathbf{M} \mathbf{w}  = 2 \mathbf{M} \mathbf{w},
\end{equation}

for all $\mathbf{M}\in \mathcal{H}_p$. For the WLS cost, the derivation is a little bit tricky and leads to

\begin{equation}
\nabla f^{\rm WLS}_{\tilde{\mathbf{\Sigma}}} = 4 {\rm diag}(\mathbf{M} {\rm diag}(\mathbf{w})^H \mathbf{M}-\mathbf{M}),
\end{equation}

where $\mathbf{M}=\tilde{\boldsymbol{\Sigma}}^{-1}{\rm diag}(\mathbf{w})\tilde{\boldsymbol{\Sigma}}$ and diag gives a diagonal matrix when the entry is a vector and a vector when the entry is a diagonal matrix.

Riemannian conjugate gradient on $\mathbb{T}_p$

[1]
Entry: Objective $f$ of problem (<ref>), starting point $\mathbf{w}_0$
Compute $\boldsymbol{\eta}_t = - {\rm grad} f(\mathbf{w}_t)$ with (<ref>)
Compute $\alpha_t$ and $\beta_t$ with [67]
Set direction $\boldsymbol{\xi}_t = \alpha_t \boldsymbol{\eta}_t + \beta_t \mathcal{T}_{\mathbf{w}_{t}\leftarrow\mathbf{w}_{t-1}}^{\mathbb{T}_p} (  \boldsymbol{\eta}_{t-1}) $
Update $\mathbf{w}_t = \mathcal{R}_{\mathbf{w}_t} ({\boldsymbol{\xi}}_t)$ with (<ref>)
Output: $\hat{\mathbf{w}}_{\boldsymbol{\theta}} = \mathbf{w}_{\rm end} \in \mathbb{T}_p$

§ LINKS WITH EXISTING ALGORITHMS

§.§ Existing algorithms as COFI-PL instances

An overview relating state-of-the-art IPL algorithms to module options of COFI-PL is presented in Table <ref>.

We can notice the prevalence of KL-IPL formulations due to their interesting Gaussian maximum likelihood roots.

In this scope, a subtlety requires some discussion:

as it is formalized in this paper, COFI-PL aims at fitting the phases of the plug-in (regularized) estimator $\tilde{\mathbf{\Sigma}}$ to its own modulus $\tilde{\mathbf{\Psi}} \overset{\Delta}{=} {\rm mod}(\tilde{\mathbf{\Sigma}})$.

For KL-IPL, this means minimizing the quadratic form on $\mathbb{T}_p$ in (<ref>), with the matrix  $\mathbf{M}_{\rm KL}= \tilde{\mathbf{\Psi}}^{-1} \circ \tilde{\mathbf{\Sigma}} $.

On the other hand, traditional formulations of KL-IPL arise from a Gaussian maximum likelihood approach.

In this case, $\tilde{{\mathbf{\Sigma}}} = \mathbf{S}$ is set by the construction of the likelihood function, and the regularization is only applied to the plug-in estimate of the unknown modulus $\tilde{\mathbf{\Psi}}$ (before inversion).

Such construction of the objective function is denoted as KL$_{\rm ML}$ in the table to highlight the distinction.

Though not studied in this paper for the sake of clarity, this difference suggests a more flexible implementation, in which one can dissociate (mix-and-match) regularization processes applied to $\tilde{\mathbf{\Psi}}^{-1} $ and $\tilde{\mathbf{\Sigma}}$ when constructing $\mathbf{M}_{\rm KL}$.

Method      Plug-in $\hat{\mathbf{\Sigma}}$     Regularization     Fitting cost     Optimization     Remark
PL [1]     SCM         KL     Elementwise on $\{ \theta_q \}_{q=1}^p$

PTA [2]     SCM         KL     BFGS on $\boldsymbol{\theta}$

CAESAR [3]      SCorr      LR (<ref>), $r=1$     KL$_{\rm ML}$

Cao et al. [4]     SCorr      LR (optional)     KL          LR uses KL$_{\rm ML}$ as in [3]

Bootstrapping [5]     Boot. Corr          KL$_{\rm ML}$          Bootstrapping method for coherence estimation

EMI  [6]     SCM          KL     EVD relax.

TMLE [68]     SCM      shrink.     (<ref>)

Zwieback [8]     SCM     shrink. and/or Tap.       KL$_{\rm ML}$

LS-PL [69]     SCM         LS      MM on $\mathbf{w}_{\boldsymbol{\theta}}$

LaMIE [10]     SCM     shrink. and Tap.      LS     BFGS on $\boldsymbol{\theta}$ or EVD relax.     Parameter selection with [54]

State-of-the-art IPL algorithms as instances of COFI-PL

§.§ Other maximum likelihood approaches

Section <ref> discussed the use of structure constrained covariance matrix estimators.

It naturally raises the following questions: is it possible to directly perform structured maximum likelihood estimation with the structure constraint in (<ref>)?

Unfortunately, the modulus-argument decomposition is not an holomorphic function [70, 71], so it is not well suited to a differentiation-based optimization process.

To overcome this issue, several works [24, 68, 25] considered the alternative decomposition $\mathbf{\Sigma} = \boldsymbol{\Psi}_{\mathcal{R}}  \circ (\mathbf{w}_{\boldsymbol{\theta}} \mathbf{w}_{\boldsymbol{\theta}}^H)$, where $\boldsymbol{\Psi}_{\mathcal{R}} \in \mathcal{S}_p^{++}$ is a real positive definite matrix.

This decomposition involves parameters in smooth manifolds, so it can be efficiently accounted for in a constrained optimization problem.

It also coincides with the modulus-argument when all entries in $\boldsymbol{\Psi}_{\mathcal{R}}$ are positive.

However, it is less restrictive, and the resulting solutions can potentially include negative entries in $\boldsymbol{\Psi}_{\mathcal{R}}$, i.e., some ambiguities that break the phase closure property (possible correction procedures are discussed in [25]).

In the Gaussian case, [68] relies on a relaxing simplification that allows to compress the likelihood in the objective function into a single determinant:

\begin{equation} \label{eq:mindetTMLE}
\begin{array}{c l}
\underset{{\mathbf{w}_{\boldsymbol{\theta}}}\in\mathbb{T}_p}{\rm minimize}
\left| \mathfrak{Re} \{ {\rm diag}(\mathbf{w}_{\boldsymbol{\theta}})^H \tilde{\mathbf{S}} {\rm diag}(\mathbf{w}_{\boldsymbol{\theta}}) \}  \right|
\end{array}
\end{equation}

Alternatively, majorization-minimization algorithms for the joint optimization of $\boldsymbol{\Psi}_{\mathcal{R}} $ and $ \mathbf{w}_{\boldsymbol{\theta}}$ are derived in [24] for the Gaussian case, and extended to rank-constrained and/or scaled Gaussian models in [25].

In any case, the term “true maximum likelihood” needs to be handled with caution because [68] appears to employ a relaxation of the Gaussian log-likelihood. Moreover, all works [24, 68, 25] rely on a matrix decomposition that is more permissive than the expected structure (<ref>) (though evidenced to be useful in practice).

§ REAL CASE STUDY: SUBSIDENCE OF MEXICO CITY

COFI-PL with KL fitting applied to Mexico InSAR dataset with various plug-in estimators (from left to right: sample covariance matrix (<ref>), Tyler's estimator (<ref>) with $u_{T}(t) = p/t$, phase-only sample covariance matrix (<ref>)) and various regularization processes (from top to bottom: no regularization, low-rank approximation (<ref>) with $k=1$, tapering (<ref>) with bandwidth $b=9$, shrinkage to identity (<ref>) with $\beta=0.1$).

COFI-PL with LS fitting applied to Mexico InSAR dataset with various plug-in estimators (from left to right: sample covariance matrix (<ref>), Tyler's estimator (<ref>) with $u_{T}(t) = p/t$, phase-only sample covariance matrix (<ref>)) and various regularization processes (from top to bottom: no regularization, low-rank approximation (<ref>) with $k=1$, tapering (<ref>) with bandwidth $b=9$, shrinkage to identity (<ref>) with $\beta=0.1$).

COFI-PL with WLS fitting applied to Mexico InSAR dataset with various plug-in estimators (from left to right: sample covariance matrix (<ref>), Tyler's estimator (<ref>) with $u_{T}(t) = p/t$, phase-only sample covariance matrix (<ref>)) and various regularization processes (from top to bottom: no regularization, low-rank approximation (<ref>) with $k=1$, tapering (<ref>) with bandwidth $b=9$, shrinkage to identity (<ref>) with $\beta=0.1$).

§.§ Dataset

To assess the performances of all proposed algorithms, we investigate one dataset with a temporal coverage of a year. The dataset consists of 31 Sentinel-1 SLC SAR images acquired between July 2019 and June 2020 over the southern region of Mexico Valley, within an endoreic basin surrounded by volcanic mountains. Pre-processing steps such as Sentinel-1 dedicated processing, coregistration, initial interferogram computation and filtering are done using the free software, SNAP, developed by ESA [72].

Results presented in the next paragraph are obtained using $p=31$ temporal samples and $n=64$ spatial samples (2D sliding window with size $8 \times 8$ pixels).

§.§ Results

For a given objective function set by the chosen matrix distance $d$, COFI-PL is applied to the data with various options of plug-in estimates and regularizations.

Fig. <ref> (resp. Fig. <ref>, and Fig. <ref>) presents the output of the algorithm when KL (resp. LS, and WLS) is used as fitting objective.

COFI-PL with KL corresponds to most of standard IPL formulations (albeit the slight difference regarding regularization discussed in Section <ref>).

In this configuration, we can clearly see in Fig. <ref> that the quality of the plug-in estimate plays an important role (cf. first line), and that the phase-only sample covariance matrix (<ref>) provides a great improvement compared to the standard sample covariance matrix.

The regularization also plays a crucial part in stabilizing the inversion of the modulus of the plug-in estimate (i.e., $\tilde{\boldsymbol{\Psi}}^{-1}$) that is required by construction with the KL fitting cost.

In this setup, the low-rank approximation (<ref>) and shrinkage to identity (<ref>) appear to be the best regularization options in order to compensate a poor estimation or a poor conditioning of this quantity.

The tapering regularization using a hard threshold appears less stable, probably because it affects the aforementioned inversion step.

Hence, other (smoother) forms of tapering could be envisioned [54].

The same conclusions can be drawn for COFI-PL with WLS fitting objective presented in Fig. <ref>.

This was to be expected because the KL and WLS distances have a similar underlying construction that involves the inverse of the modulus of the plug-in estimate.

We can still notice some differences: WLS is even more sensitive to the inversion instabilities induced by the tapering (since it depends quadraticly on $\tilde{\boldsymbol{\Psi}}^{-1}$), however it provides an interesting alternative to KL is some setups (notably when using Tyler's estimator as plug-in).

The quality and robustness of the plug-in estimate also plays a role for COFI-PL with LS fitting objective, as we can see an improvement brought by the phase-only sample covariance matrix in the first line of Fig. <ref>.

Interestingly, the LS setup appears much less sensitive to the conditioning of the plug-in estimate because no matrix inversion is required.

Hence results appear more stable with respect to the various regularization options.

A closer inspection of the outputs further reveals that regularization is still beneficial, notably when tapering is used.

As a matter of fact, LS fitting on a banded phase-only sample covariance matrix allowed to obtain the best visual results on this dataset.

In terms of optimization, majorization-minimization and Riemannian gradient descent yielded the same outputs for KL and LS fitting objectives, and WLS could only be evaluated with the Riemannian optimization framework.

In practice, majorization-minimization for KL and LS was experienced to be easier to tune (as no side parameters such as step size are involved) and provided solutions more quickly than a plain Riemannian gradient descent with constant step size.

Note that faster Riemannian optimization schemes could probably be obtained with more involved adaptive step size selection rules or alternate Riemannian optimization methods (e.g., Riemannian conjugate gradient or BFGS).

However, this last point goes beyond the scope of this paper.

§ CONCLUSION

This paper presented COFI-PL: a compact framework to design IPL algorithms based on the covariance fitting approach.

The formulation encompasses existing methods initially derived as approximate maximum likelihood, and allowed us to present some promising alternatives using various robust plug-in estimators, regularization methods, and fitting objective functions.

It also provides a clear and systematic framework to perform ablation studies, investigate the interest of alternates options in a single module, or assess the coupling effects when considering a joint estimation-regularization method of the covariance matrix.

Two methods to deal with the resulting optimization problems were then introduced: majorization-minimization and Riemannian optimization on the torus of phase-only complex vector.

Generic principles and tools were presented so that any cost optimization problem constructed within the COFI-PL framework can be addressed efficiently.

A real-world case study regarding the subsidience of Mexico city finally illustrated the interest of the proposed approach.

Experiments on this dataset notably evidenced the practical interest of using $i$) a LS fitting cost rather than the traditional KL divergence emanating from the maximum likelihood perspective; $ii$) the phase-only sample covariance matrix (a newly introduced covariance matrix estimate in the context of IPL) as plug-in estimate.

[1]
Andrea Monti Guarnieri and Stefano Tebaldini,
“On the Exploitation of Target Statistics for SAR Interferometry
IEEE Transactions on Geoscience and Remote Sensing, vol. 46,
no. 11, pp. 3436–3443, 2008.

[2]
Alessandro Ferretti, Alfio Fumagalli, Fabrizio Novali, Claudio Prati, Fabio
Rocca, and Alessio Rucci,
“A New Algorithm for Processing Interferometric Data-Stacks:
IEEE Transactions on Geoscience and Remote Sensing, vol. 49,
no. 9, pp. 3460–3470, 2011.

[3]
Gianfranco Fornaro, Simona Verde, Diego Reale, and Antonio Pauciullo,
“CAESAR: An Approach Based on Covariance Matrix Decomposition to
Improve Multibaseline–Multitemporal Interferometric SAR Processing,”
IEEE Transactions on Geoscience and Remote Sensing, vol. 53,
no. 4, pp. 2050–2065, 2015.

[4]
Ning Cao, Hyongki Lee, and Hahn Chul Jung,
“Mathematical Framework for Phase-Triangulation Algorithms in
Distributed-Scatterer Interferometry,”
IEEE Geoscience and Remote Sensing Letters, vol. 12, no. 9, pp.
1838–1842, 2015.

[5]
Mi Jiang and Andrea Monti Guarnieri,
“Distributed scatterer interferometry with the refinement of
spatiotemporal coherence,”
IEEE Transactions on Geoscience and Remote Sensing, vol. 58,
no. 6, pp. 3977–3987, 2020.

[6]
Homa Ansari, Francesco De Zan, and Richard Bamler,
“Efficient phase estimation for interferogram stacks,”
IEEE Transactions on Geoscience and Remote Sensing, vol. 56,
no. 7, pp. 4109–4125, 2018.

[7]
Dinh Ho Tong Minh and Yen-Nhi Ngo,
“Compressed SAR Interferometry in the Big Data Era,”
Remote Sensing, vol. 14, no. 2, pp. 390, 2022.

[8]
Simon Zwieback,
“Cheap, valid regularizers for improved interferometric phase
IEEE Geoscience and Remote Sensing Letters, vol. 19, pp. 1–4,

[9]
Phan Viet Hoa Vu, Arnaud Breloy, Frédéric Brigui, Yajing Yan, and
Guillaume Ginolhac,
“Covariance fitting based InSAR Phase Linking,”
in IGARSS 2023-2023 IEEE International Geoscience and Remote
Sensing Symposium. IEEE, 2023, pp. 8234–8237.

[10]
Yusong Bai, Jian Kang, Xiang Ding, Anping Zhang, Zhe Zhang, and Naoto Yokoya,
“LaMIE: Large-Dimensional Multipass InSAR Phase Estimation for
Distributed Scatterers,”
IEEE Transactions on Geoscience and Remote Sensing, vol. 61,
pp. 1–15, 2023.

[11]
Michael Schmitt, Johannes L Schönberger, and Uwe Stilla,
“Adaptive covariance matrix estimation for multi-baseline InSAR
data stacks,”
IEEE Transactions on Geoscience and Remote Sensing, vol. 52,
no. 11, pp. 6807–6817, 2014.

[12]
Yuanyuan Wang and Xiao Xiang Zhu,
“Robust estimators for multipass SAR interferometry,”
IEEE Transactions on Geoscience and Remote Sensing, vol. 54,
no. 2, pp. 968–980, 2015.

[13]
Dinh Ho Tong Minh and Stefano Tebaldini,
“Interferometric Phase Linking: Algorithm, application, and
IEEE Geoscience and Remote Sensing Magazine, vol. 11, no. 3,
pp. 46–62, 2023.

[14]
Bjorn Ottersten, Petre Stoica, and Richard Roy,
“Covariance Matching Estimation Techniques for Array Signal
Processing Applications,”
Digital Signal Processing, vol. 8, no. 3, pp. 185–210, 1998.

[15]
Shahram Shahbazpanahi, Shahrokh Valaee, and Alex B. Gershman,
“A covariance fitting approach to parametric localization of
multiple incoherently distributed sources,”
IEEE Transactions on Signal Processing, vol. 52, no. 3, pp.
592–600, 2004.

[16]
Tarik Yardibi, Jian Li, Petre Stoica, Nikolas S. Zawodny, and III Cattafesta,
Louis N.,
“A covariance fitting approach for correlated acoustic source
The Journal of the Acoustical Society of America, vol. 127, no.
5, pp. 2920–2931, 05 2010.

[17]
Petre Stoica, Prabhu Babu, and Jian Li,
“SPICE: A Sparse Covariance-Based Estimation Method for Array
IEEE Transactions on Signal Processing, vol. 59, no. 2, pp.
629–638, 2011.

[18]
Bruno Mériaux, Chengfang Ren, Mohammed Nabil El Korso, Arnaud Breloy, and
Philippe Forster,
“Robust-COMET for covariance estimation in convex structures:
Algorithm and statistical properties,”
in 2017 IEEE 7th International Workshop on Computational
Advances in Multi-Sensor Adaptive Processing (CAMSAP), 2017, pp. 1–5.

[19]
Bruno Mériaux, Chengfang Ren, Mohammed Nabil El Korso, Arnaud Breloy, and
Philippe Forster,
“Robust estimation of structured scatter matrices in (mis)matched
Signal Processing, vol. 165, pp. 163–174, 2019.

[20]
Bruno Mériaux, Chengfang Ren, Arnaud Breloy, Mohammed Nabil El Korso, Philippe
Forste, and Jean-Philippe Ovarlez,
“On the Recursions of Robust COMET Algorithm for Convexly
Structured Shape Matrix,”
in 2019 27th European Signal Processing Conference (EUSIPCO),
2019, pp. 1–5.

[21]
Bruno Mériaux, Chengfang Ren, Arnaud Breloy, Mohammed Nabil El Korso, and
Philippe Forster,
“Matched and Mismatched Estimation of Kronecker Product of Linearly
Structured Scatter Matrices Under Elliptical Distributions,”
IEEE Transactions on Signal Processing, vol. 69, pp. 603–616,

[22]
David R Hunter and Kenneth Lange,
“A tutorial on MM algorithms,”
The American Statistician, vol. 58, no. 1, pp. 30–37, 2004.

[23]
Ying Sun, Prabhu Babu, and Daniel P Palomar,
“Majorization-minimization algorithms in signal processing,
communications, and machine learning,”
IEEE Transactions on Signal Processing, vol. 65, no. 3, pp.
794–816, 2016.

[24]
Phan Viet Hoa Vu, Frédéric Brigui, Arnaud Breloy, Yajing Yan, and Guillaume
“A New Phase Linking Algorithm for Multi-temporal InSAR based on
the Maximum Likelihood Estimator,”
in IGARSS 2022 - 2022 IEEE International Geoscience and Remote
Sensing Symposium, 2022, pp. 76–79.

[25]
Phan Viet Hoa Vu, Arnaud Breloy, Frédéric Brigui, Yajing Yan, and Guillaume
“Robust Phase Linking in InSAR,”
IEEE Transactions on Geoscience and Remote Sensing, vol. 61,
pp. 1–11, 2023.

[26]
Steven T Smith,
“Optimum phase-only adaptive nulling,”
IEEE transactions on signal processing, vol. 47, no. 7, pp.
1835–1843, 1999.

[27]
Mohamed A ElMossallamy, Karim G Seddik, Wei Chen, Li Wang, Geoffery Ye Li, and
Zhu Han,
“Ris optimization on the complex circle manifold for interference
mitigation in interference channels,”
IEEE Transactions on Vehicular Technology, vol. 70, no. 6, pp.
6184–6189, 2021.

[28]
Weijie Xiong, Jinfeng Hu, and Kai Zhong,
“Mimo radar transmit beampattern matching via manifold
in ICASSP 2023-2023 IEEE International Conference on Acoustics,
Speech and Signal Processing (ICASSP). IEEE, 2023, pp. 1–5.

[29]
P-A Absil, Robert Mahony, and Rodolphe Sepulchre,
Optimization algorithms on matrix manifolds,
Princeton University Press, 2008.

[30]
Nicolas Boumal,
An introduction to optimization on smooth manifolds,
Cambridge University Press, 2023.

[31]
Xinping Deng, Carlos López-Martínez, Jinsong Chen, and Pengpeng Han,
“Statistical modeling of polarimetric SAR data: A survey and
Remote Sensing, vol. 9, no. 4, pp. 348, 2017.

[32]
Ammar Mian, Guillaume Ginolhac, Jean-Philippe Ovarlez, and
Abdourrahmane Mahamane Atto,
“New robust statistics for change detection in time series of
multivariate SAR images,”
IEEE Transactions on Signal Processing, vol. 67, no. 2, pp.
520–534, 2018.

[33]
Esa Ollila, David E Tyler, Visa Koivunen, and H Vincent Poor,
“Complex elliptically symmetric distributions: Survey, new results
and applications,”
IEEE Transactions on signal processing, vol. 60, no. 11, pp.
5597–5625, 2012.

[34]
Ricardo Antonio Maronna,
“Robust M-estimators of multivariate location and scatter,”
The annals of statistics, pp. 51–67, 1976.

[35]
David E. Tyler,
“A Distribution-Free $M$-Estimator of Multivariate Scatter,”
The Annals of Statistics, vol. 15, no. 1, pp. 234 – 251, 1987.

[36]
Georgy Shevlyakov and Pavel Smirnov,
“Robust estimation of the correlation coefficient: An attempt of
Austrian Journal of Statistics, vol. 40, no. 1&2, pp.
147–156, 2011.

[37]
Georgy L Shevlyakov and Hannu Oja,
Robust correlation: Theory and applications, vol. 3,
John Wiley & Sons, 2016.

[38]
Yonatan Woodbridge, Gal Elidan, and Ami Wiesel,
“Signal detection in complex structured para normal noise,”
IEEE Transactions on Signal Processing, vol. 65, no. 9, pp.
2306–2316, 2017.

[39]
Olivier Ledoit and Michael Wolf,
“A well-conditioned estimator for large-dimensional covariance
Journal of Multivariate Analysis, vol. 88, no. 2, pp. 365–411,

[40]
Yilun Chen, Ami Wiesel, and Alfred O. Hero,
“Robust Shrinkage Estimation of High-Dimensional Covariance
IEEE Transactions on Signal Processing, vol. 59, no. 9, pp.
4097–4107, 2011.

[41]
Esa Ollila and Elias Raninen,
“Optimal shrinkage covariance matrix estimation under random
sampling from elliptical distributions,”
IEEE Transactions on Signal Processing, vol. 67, no. 10, pp.
2707–2719, 2019.

[42]
Esa Ollila, Daniel P Palomar, and Frédéric Pascal,
“Shrinking the eigenvalues of M-estimators of covariance matrix,”
IEEE Transactions on Signal Processing, vol. 69, pp. 256–269,

[43]
Ammar Mian, Antoine Collas, Arnaud Breloy, Guillaume Ginolhac, and
Jean-Philippe Ovarlez,
“Robust low-rank change detection for multivariate SAR image time
IEEE Journal of Selected Topics in Applied Earth Observations
and Remote Sensing, vol. 13, pp. 3545–3556, 2020.

[44]
Petre Stoica and Yngve Selen,
“Model-order selection: a review of information criterion rules,”
IEEE Signal Processing Magazine, vol. 21, no. 4, pp. 36–47,

[45]
Ian T Jolliffe,
“Principal component analysis,”
Technometrics, vol. 45, no. 3, pp. 276, 2003.

[46]
Michael E Tipping and Christopher M Bishop,
“Probabilistic principal component analysis,”
Journal of the Royal Statistical Society. Series B (Statistical
Methodology), vol. 61, no. 3, pp. 611–622, 1999.

[47]
Bosung Kang, Vishal Monga, and Muralidhar Rangaswamy,
“Rank-constrained maximum likelihood estimation of structured
covariance matrices,”
IEEE Transactions on Aerospace and Electronic Systems, vol. 50,
no. 1, pp. 501–515, 2014.

[48]
Homa Ansari, Francesco De Zan, and Alessandro Parizzi,
“Study of systematic bias in measuring surface deformation with SAR
IEEE Transactions on Geoscience and Remote Sensing, vol. 59,
no. 2, pp. 1285–1301, 2020.

[49]
Peng Shen, Changcheng Wang, and Bei An,
“AdpPL: An adaptive phase linking-based distributed scatterer
interferometry with emphasis on interferometric pair selection optimization
and adaptive regularization,”
Remote Sensing of Environment, vol. 295, pp. 113687, 2023.

[50]
T Tony Cai, Cun-Hui Zhang, and Harrison H Zhou,
“Optimal rates of convergence for covariance matrix estimation,”
The Annals of Statistics, vol. 38, no. 4, pp. 2118 – 2144,

[51]
Joseph R Guerci,
“Theory and application of covariance matrix tapers for robust
adaptive beamforming,”
IEEE Transactions on Signal Processing, vol. 47, no. 4, pp.
977–985, 1999.

[52]
Peter J Bickel and Elizaveta Levina,
“Regularized estimation of large covariance matrices,”
The Annals of Statistics, vol. 36, no. 1, pp. 199 – 227, 2008.

[53]
Peter J Bickel and Elizaveta Levina,
“Covariance regularization by thresholding,”
The Annals of Statistics, vol. 36, no. 6, pp. 2577 – 2604,

[54]
Esa Ollila and Arnaud Breloy,
“Regularized tapered sample covariance matrix,”
IEEE Transactions on Signal Processing, vol. 70, pp.
2306–2320, 2022.

[55]
Ami Wiesel,
“Unified framework to regularized covariance estimation in scaled
Gaussian models,”
IEEE Transactions on Signal Processing, vol. 60, no. 1, pp.
29–38, 2011.

[56]
Ying Sun, Prabhu Babu, and Daniel P Palomar,
“Regularized Tyler's scatter estimator: Existence, uniqueness, and
IEEE Transactions on Signal Processing, vol. 62, no. 19, pp.
5143–5156, 2014.

[57]
Frédéric Pascal, Yacine Chitour, and Yihui Quek,
“Generalized robust shrinkage estimator and its application to
STAP detection problem,”
IEEE Transactions on Signal Processing, vol. 62, no. 21, pp.
5640–5651, 2014.

[58]
Esa Ollila and David E Tyler,
“Regularized $ M $-estimators of scatter matrix,”
IEEE Transactions on Signal Processing, vol. 62, no. 22, pp.
6059–6070, 2014.

[59]
Markus Even and Karsten Schulz,
“InSAR deformation analysis with distributed scatterers: A review
complemented by new advances,”
Remote Sensing, vol. 10, no. 5, pp. 744, 2018.

[60]
Ying Sun, Prabhu Babu, and Daniel P Palomar,
“Robust estimation of structured covariance matrix for heavy-tailed
elliptical distributions,”
IEEE Transactions on Signal Processing, vol. 64, no. 14, pp.
3576–3590, 2016.

[61]
Zhang Yunjun, Heresh Fattahi, and Falk Amelung,
“Small baseline InSAR time series analysis: Unwrapping error
correction and noise reduction,”
Computers & Geosciences, vol. 133, pp. 104331, 2019.

[62]
A Monti Guarnieri and Stefano Tebaldini,
“Hybrid Cramér–Rao bounds for crustal displacement field
estimators in SAR interferometry,”
IEEE signal processing letters, vol. 14, no. 12, pp.
1012–1015, 2007.

[63]
Jiyuan Hu, Wenhao Wu, Mahdi Motagh, Fen Qin, Jiayao Wang, Shangyi Pan, Jiming
Guo, and Chunyu Zhang,
“FIM-based DSInSAR method for mapping and monitoring of reservoir
bank landslides: an application along the Lancang River in China,”
Landslides, vol. 20, no. 11, pp. 2479–2495, 2023.

[64]
Meisam Razaviyayn, Mingyi Hong, and Zhi-Quan Luo,
“A Unified Convergence Analysis of Block Successive Minimization
Methods for Nonsmooth Optimization,”
SIAM on Optimization, vol. 23, no. 2, pp. 1126–1153, 2013.

[65]
Mojtaba Soltanalian and Petre Stoica,
“Designing unimodular codes via quadratic optimization,”
IEEE Transactions on Signal Processing, vol. 62, no. 5, pp.
1221–1234, 2014.

[66]
Arnaud Breloy, Sandeep Kumar, Ying Sun, and Daniel P Palomar,
“Majorization-minimization on the Stiefel manifold with application
to robust sparse PCA,”
IEEE Transactions on Signal Processing, vol. 69, pp.
1507–1520, 2021.

[67]
Magnus R Hestenes and Eduard Stiefel,
“Methods of conjugate gradients for solving linear systems,”
Journal of research of the National Bureau of Standards, vol.
49, no. 6, pp. 409–436, 1952.

[68]
Chisheng Wang, Xiang-Sheng Wang, Yaping Xu, Bochen Zhang, Mi Jiang, Siting
Xiong, Qin Zhang, Weidong Li, and Qingquan Li,
“A New Likelihood Function for Consistent Phase Series Estimation
in Distributed Scatterer Interferometry,”
IEEE Transactions on Geoscience and Remote Sensing, vol. 60,
pp. 1–14, 2022.

[69]
Phan Viet Hoa Vu, Arnaud Breloy, Frédéric Brigui, Yajing Yan, and
Guillaume Ginolhac,
“Covariance fitting based InSAR Phase Linking,”
in IGARSS 2023 - 2023 IEEE International Geoscience and Remote
Sensing Symposium, 2023.

[70]
Walter Rudin,
Real and Complex Analysis, 3rd Ed.,
McGraw-Hill, Inc., 1987.

[71]
Walter Appel and Emmanuel Kowalski,
Mathematics for physics and physicists, vol. 47,
Princeton University Press Princeton, NJ, 2007.

[72]
ESA (European Space Agency),
“SNAP,” <https://step.esa.int/main/toolboxes/snap>.