Skip to content
BY 4.0 license Open Access Published by De Gruyter January 5, 2022

Diffusion tensor regularization with metric double integrals

Leon Frischauf ORCID logo, Melanie Melching ORCID logo and Otmar Scherzer ORCID logo

Abstract

In this paper, we propose a variational regularization method for denoising and inpainting of diffusion tensor magnetic resonance images. We consider these images as manifold-valued Sobolev functions, i.e. in an infinite dimensional setting, which are defined appropriately. The regularization functionals are defined as double integrals, which are equivalent to Sobolev semi-norms in the Euclidean setting. We extend the analysis of [14] concerning stability and convergence of the variational regularization methods by a uniqueness result, apply them to diffusion tensor processing, and validate our model in numerical examples with synthetic and real data.

MSC 2010: 47A52; 65J20

1 Introduction

In this paper, we investigate denoising and inpainting of diffusion tensor (magnetic resonance) images (DTMRI) with a derivative-free, non-local variational regularization technique proposed, implemented and analyzed first in [14].

The proposed regularization functionals generalize equivalent definitions of the Sobolev semi-norms, which have been derived in the fundamental work of [11] and follow-up work [16, 40]. These papers provide a derivative-free representation of Sobolev semi-norms for intensity and vector-valued functions. The beauty of this representation is that it allows for a straightforward definition of Sobolev energies of manifold-valued data (see [14]), without requiring difficult differential geometric concepts (as for instance in [21, 15]).

Diffusion tensor images are considered to be re-presentable as functions from an image domain Ω R n , with n = 2 , 3 , respectively, into the manifold of symmetric, positive definite matrices in R m × m , denoted by 𝐾 in the following – for DTMRI images m = 3 . Therefore, they are ideal objects to check the efficiency of the proposed regularization techniques. A measured diffusion tensor is often very noisy, and post-processing steps for noise removal are important. Even more, due to the noise, it is possible that a measured tensor has negative eigenvalues, which is not physical, and thus often the whole tensor at this point is omitted, leading to incomplete data. Then the missing information has to be inpainted before visualization.

Variational regularization of vector and matrix manifold-valued functions has been considered before, for instance in [49, 44, 46, 5, 8, 30, 50, 36] and [9, 12, 37, 34]. Non-local regularization formulations were studied for example in [10, 28] for filtering tensor-valued functions, and also [22, 29] for filtering intensity images. An overview of diffusion and regularization techniques for vector- and matrix-valued data is given in [49].

Variational methods for denoising and inpainting attempt to find a good compromise between matching some given noisy, tensor-valued data w δ : Ω K and prior information on the desired solution w 0 : Ω K , also called noise-free or ideal solution. The choice of prior knowledge on w 0 is

  1. that it is an element of the set W s , p ( Ω ; K ) , with a metric 𝑑 on 𝐾, the set of positive definite matrices, which is a subset of the fractional Sobolev space W s , p ( Ω ; R m × m ) , with s ( 0 , 1 ) and p ( 1 , ) ,

  2. and that

    (1.1) Φ [ d ] l ( w ) := Ω × Ω d p ( w ( x ) , w ( y ) ) | x - y | n + p s ρ l ( x - y ) d ( x , y )

    is relatively small. The function 𝜌 is a non-negative and radially symmetric mollifier with an on-off indicator l { 0 , 1 } denoting whether the mollifier is used or not. Note that, in case that d = d R m × m is the Euclidean metric and if we choose in addition l = 0 , Φ [ d R m × m ] 0 becomes the fractional Sobolev semi-norm.

The compromise of approximating w δ with a function in W s , p ( Ω ; K ) with a small energy term Φ [ d ] l ( w ) is achieved by minimization of the functional

(1.2) F [ d ] α , w δ ( w ) := Ω χ Ω D ( x ) d p ( w ( x ) , w δ ( x ) ) d x + α Φ [ d ] l ( w ) ,

where the parameter α > 0 determines the preference of staying close to the given function w δ in Ω D and a small energy Φ [ d ] l ( w ) . One should not confuse the energy term with double integral representations approximating semi-norms on manifolds (see for instance [24, 27, 19]).

The indicator function of Ω D ,

χ Ω D ( x ) = { 1 if x Ω D , 0 otherwise ,

used in (1.2) allows us to consider the two tasks of denoising ( D = ) and inpainting ( D ) within one analysis.

The paper is organized as follows. In Section 2, we constitute our notation and setting used to analyze variational methods for DTMRI data processing. We review regularization results from [14] in Section 3. In Section 4, we verify that these results from Section 3 are applicable in the context of diffusion tensor imaging, meaning that we show that the functional F [ d ] α , w δ defined in (1.2) attains a minimizer and fulfills a stability as well as a convergence result. Furthermore, we extend the analysis and give a uniqueness result using differential geometric properties of symmetric, positive definite matrices, where it is of particular importance that these matrices endowed with the log-Euclidean metric form a flat Hadamard manifold. In Section 5, we give more details on the numerical minimization of the regularization functional, and discuss different variants. In the last section, Section 6, we show numerical results for denoising and inpainting problems of synthetic and real DTMRI data.

2 Notation and setting

In the beginning, we summarize basic notation and assumptions used throughout the paper. In the theoretical part, we work with general dimensions n , m N , while we consider the particular case n = 2 , m = 3 , that is two-dimensional slices of a three-dimensional DTMRI image, in the numerical examples in Section 6.

Assumption 2.1

We assume the following.

  1. Ω R n is a nonempty, bounded and connected open set with Lipschitz boundary, and D Ω is measurable.

  2. p ( 1 , ) , s ( 0 , 1 ) and l { 0 , 1 } .

  3. K R m × m is a nonempty and closed subset of R m × m .

  4. d R m × m : R m × m × R m × m [ 0 , ) denotes the Euclidean distance induced by the (Frobenius norm) on R m × m .

  5. d := d K : K × K [ 0 , ) denotes an arbitrary metric on 𝐾 which is equivalent to d R m × m .

Moreover, we need the definition of a mollifier which appears in the regularizer of the functional in (1.2).

Definition 2.2

Definition 2.2 (Mollifier)

We call ρ C c ( R n ; R ) a mollifier if

  • 𝜌 is a non-negative, radially symmetric function,

  • R n ρ ( x ) d x = 1 and

  • there exists some 0 < τ < ρ L ( R n ; R ) and η := η τ > 0 such that { z R n : ρ ( z ) τ } = { z R n : | z | η } .

The last condition holds for instance if 𝜌 is radially decreasing satisfying ρ ( 0 ) > 0 .

2.1 Diffusion tensors

It is commonly assumed that the recorded diffusion tensor images are functions with values which are symmetric, positive definite matrices. Hence we make the assumption that

w , w δ : Ω SPD ( m ) ,

where SPD ( m ) is the set of symmetric, positive definite, real m × m matrices defined below in (2.1). When working with data from MRI measurements, m = 3 .

In the following definition, we summarize sets of matrices and associated norms on the sets.

Definition 2.3

The vector space of symmetric matrices is

SYM ( m ) := { M R m × m : M T = M } .

Additionally, we define the set of symmetric, positive definite m × m matrices

(2.1) SPD ( m ) := { M SYM ( m ) : x T M x > 0 for x R m { 0 } } .

The set of symmetric, positive definite matrices with bounded spectrum is

(2.2) SPD [ ε ¯ , ε ¯ ] spec ( m ) := { M SPD ( m ) : spec ( M ) [ ε ¯ , ε ¯ ] } ,

where spec denotes the spectrum of a given matrix. For diffusion tensors, the spectrum is real.

The set of symmetric, positive definite matrices with bounded logarithm is

(2.3) SPD z Log ( m ) := { M SPD ( m ) : Log ( M ) F z } ,

where Log is the matrix logarithm defined later in Definition 4.2(ii) and F denotes the Frobenius norm defined as

(2.4) M F = i , j = 1 m | m i j | 2 .

When working with DTMRI data, in particular in Section 6, we will chose K = SPD z Log ( 3 ) . In the general theory stated in Section 3, any nonempty and bounded set can be taken.

From now on and whenever possible, we omit the space dimension and write SYM , SPD , SPD [ ε ¯ , ε ¯ ] spec and SPD z Log instead of SYM ( m ) , SPD ( m ) , SPD [ ε ¯ , ε ¯ ] spec ( m ) and SPD z Log ( m ) .

2.2 Fractional Sobolev spaces

Moreover, we need the definition of fractional Sobolev spaces and associated subsets.

Definition 2.4

Definition 2.4 (Sobolev spaces of fractional order)

Let Assumption 2.1 hold.

  • We denote by L p ( Ω ; R m × m ) the Lebesgue space of matrix-valued functions.

  • The Sobolev space W 1 , p ( Ω ; R m × m ) consists of all weakly differentiable functions in L p ( Ω ; R m × m ) for which

    w W 1 , p ( Ω ; R m × m ) := ( w L p ( Ω ; R m × m ) p + Ω w ( x ) F p d x ) 1 / p < ,

    where w is the Jacobian of 𝑤 and | w | W 1 , p ( Ω ; R m × m ) := ( Ω w ( x ) F p d x ) 1 / p is the Sobolev semi-norm.

  • The fractional Sobolev space of order 𝑠 is defined (cf. [1]) as the set

    W s , p ( Ω ; R m × m ) := { w L p ( Ω ; R m × m ) : w ( x ) - w ( y ) F | x - y | n / p + s L p ( Ω × Ω ; R ) }

    equipped with the norm

    w W s , p ( Ω ; R m × m ) := ( w L p ( Ω ; R m × m ) p + | w | W s , p ( Ω ; R m × m ) p ) 1 / p ,

    where | w | W s , p ( Ω ; R m × m ) is the semi-norm on W s , p ( Ω ; R m × m ) , defined by

    | w | W s , p ( Ω ; R m × m ) := ( Ω × Ω w ( x ) - w ( y ) F p | x - y | n + p s d ( x , y ) ) 1 / p for all w W s , p ( Ω ; R m × m ) .

  • We define the fractional Sobolev set of order 𝑠 with data in 𝐾 as

    (2.5) W s , p ( Ω ; K ) := { w W s , p ( Ω ; R m × m ) : w ( x ) K for a.e. x Ω } .

    The Lebesgue set with data in 𝐾 is defined as

    L p ( Ω ; K ) := { w L p ( Ω ; R m × m ) : w ( x ) K for a.e. x Ω } .

Note that L p ( Ω ; K ) and W s , p ( Ω ; K ) are sets and not linear spaces because summation of elements in 𝐾 is typically not closed in 𝐾.

3 Metric double integral regularization on closed subsets

We start this section by stating conditions under which the regularization functional in (1.2) attains a minimizer and fulfills a stability as well as a convergence result. Therefore, we recall results established in [14]. There the authors define a regularization functional inspired by the work of Bourgain, Brézis and Mironescu [11, 40, 16]. The analysis in turn is based on [41]. We apply these results to diffusion tensor image denoising and inpainting in the next section.

We start by stating general conditions on the exact data w 0 , the noisy data w δ and the functional F [ d ] α , w δ , defined in (1.2).

Assumption 3.1

Let Assumption 2.1 hold. Moreover, let w 0 , w δ L p ( Ω ; K ) , and let 𝜌 be a mollifier as defined in Definition 2.2. We assume that

  1. for every t > 0 and α > 0 , the level sets

    level ( F [ d ] α , w 0 ; t ) := { w W s , p ( Ω ; K ) : F [ d ] α , w 0 ( w ) t }

    are weakly sequentially pre-compact in W s , p ( Ω ; R m × m ) .

  2. There exists t ¯ > 0 such that level ( F [ d ] α , w 0 ; t ¯ ) is nonempty.

Remark 3.2

If Assumption 2.1 is fulfilled and in particular when performing image denoising ( D = ) or inpainting ( D ) of functions with values in 𝐾, then the functional (1.2) with Φ [ d ] l as in (1.1) defined on W s , p ( Ω ; K ) satisfies Assumption 3.1 (cf. [14]).

According to [14], we now have the following result giving existence of a minimizer of the functional in (1.2) as well as a stability and convergence result.

Theorem 3.3

Let Assumption 3.1 hold (which is guaranteed by Remark 3.2). For the functional defined in (1.2) over W s , p ( Ω ; K ) with Φ [ d ] l defined in (1.1), the following results hold.

  • Existence: For every v L p ( Ω ; K ) and α > 0 , the functional F [ d ] α , v : W s , p ( Ω ; K ) [ 0 , ) attains a minimizer in W s , p ( Ω ; K ) .

  • Stability: Let α > 0 be fixed, w δ L p ( Ω ; K ) , and let ( v k ) k N be a sequence in L p ( Ω ; K ) such that

    w δ - v k L p ( Ω ; K ) 0 .

    Then every sequence ( w k ) k N satisfying

    w k arg min { F [ d ] α , v k ( w ) : w W s , p ( Ω ; K ) }

    has a converging subsequence with respect to the weak topology of W s , p ( Ω ; R m × m ) . The limit w ~ of every such converging subsequence ( w k j ) j N is a minimizer of F [ d ] α , w δ . Moreover, ( Φ [ d ] l ( w k j ) ) j N converges to Φ [ d ] l ( w ~ ) .

  • Convergence: Let α : ( 0 , ) ( 0 , ) be a function satisfying α ( δ ) 0 and δ p α ( δ ) 0 for δ 0 . Let ( δ k ) k N be a sequence of positive real numbers converging to 0. Moreover, let ( v k ) k N be a sequence in L p ( Ω ; K ) with w 0 - v k L p ( Ω ; K ) δ k , and set α k := α ( δ k ) . Then every sequence ( w k ) k N defined as

    w k arg min { F [ d ] α k , v k ( w ) : w W s , p ( Ω ; K ) }

    has a weakly converging subsequence w k j w 0 as j (with respect to the topology of W s , p ( Ω ; R m × m ) ). In addition, Φ [ d ] l ( w k j ) Φ [ d ] l ( w 0 ) . Moreover, it follows that even w k w 0 weakly (with respect to the topology of W s , p ( Ω ; R m × m ) ) and Φ [ d ] l ( w k ) Φ [ d ] l ( w 0 ) .

In the theorem above, stability (with respect to the L p -norm) ensures that the minimizers of F [ d ] α , w δ depend continuously on the given data w δ . We emphasize that, in an Euclidean setting (that is on W s , p ( Ω , R m × m ) , for s > 0 and p > 1 , one could expect convergence in even stronger norms. However, here, we have to make sure that the traces into K R m × m are well-defined in appropriate Sobolev spaces, which requires additional compactness assumptions, or in other words, stronger regularization.

In the next section, we apply Theorem 3.3 to diffusion tensor images, i.e. when choosing 𝐾 as a closed subset of the symmetric, positive definite matrices.

4 Diffusion tensor regularization

The goal of this section is to define appropriate fractional order Sobolev sets as defined in (2.5) of functions which can represent diffusion tensor images. To this end, we use the set of symmetric, positive definite m × m matrices with bounded logarithm (defined in (2.3))

(4.1) K = SPD z Log

and associate it with the log-Euclidean metric, defined below in (4.6). This metric was shown to be an adequate distance measure for DTMRI, see e.g. [20, 2].

Below, we show that Theorem 3.3 applies to the regularization functional in (1.2) with the particular choice K = SPD z Log . In addition to what follows from the general theory from [14] in a straightforward manner, we present a uniqueness result for the minimizer of the regularization functional.

We begin by defining needed concepts from matrix calculus. When working with symmetric, positive definite matrices, many of the operations below reduce to their scalar counterpart applied to the eigenvalues.

4.1 Matrix calculus

We start this section by repeating basic definitions known from matrix calculus (see for instance [35]). Especially the matrix logarithm is needed to define the log-Euclidean metric on the symmetric, positive definite matrices.

Lemma 4.1

Lemma 4.1 (Matrix properties)

The following statements hold.

  1. Eigendecomposition: Let A SYM with eigenvalues ( λ i ) i = 1 m . Then we can write

    A = U Λ U T ,

    where U R m × m is the orthonormal matrix whose 𝑖-th column consists of the 𝑖-th normalized eigenvector of 𝐴. Hence we have that U U T = 1 m , where 1 m denotes the identity matrix in R m × m . Further, Λ is the diagonal matrix whose diagonal entries are the corresponding eigenvalues, Λ = diag ( λ 1 , , λ m ) .

  2. If U , V R m × m are both unitary, then so are U V T , U T V , V U T and V T U .

Next we state the definitions of the matrix exponential and logarithm, see in particular [38, 3].

Definition 4.2

Let A , B SYM with corresponding eigendecompositions A = U Λ A U T and B = V Λ B V T , where U , V R m × m unitary and Λ A = diag ( λ 1 , , λ m ) , Λ B = diag ( μ 1 , , μ m ) R m × m diagonal.

  1. Exponential map: The exponential map is defined as

    Exp ( A ) = Exp ( U Λ A U T ) = U Exp ( Λ A ) U T .

    It holds that

    Exp ( Λ A ) = diag ( e λ 1 , , e λ m ) ,

    where e : R R 0 denotes the (scalar) exponential function, and Exp : SYM SPD is a diffeomorphism [3, Theorem 2.8].

  2. Logarithm: If Exp ( A ) = B , then A is the matrix logarithm of 𝐵. It is defined as

    Log ( B ) = Log ( V Λ B V T ) = V Log ( Λ B ) V T .

    Moreover,

    Log ( Λ B ) = diag ( log ( μ 1 ) , , log ( μ m ) ) ,

    where log : R 0 R is the (scalar) natural logarithm, i.e. log := log e .

    When restricting to symmetric, positive definite matrices, Log : SPD SYM is a diffeomorphism [3, Theorem 2.8].

The previous definition, Definition 4.2, shows that the exponential and logarithm of a symmetric (positive definite) matrix can be computed easily due to the eigendecomposition (see Remark 4.1) by calculating the scalar exponential map and logarithm of the eigenvalues.

Remark 4.3

Remark 4.3 (Matrix logarithm)

For a general matrix in R m × m , the matrix logarithm is not unique. Matrices with positive eigenvalues have a unique, symmetric logarithm, called the principal logarithm [3].

The next lemma states properties of the Frobenius norm (recall (2.4)).

Lemma 4.4

Lemma 4.4 (Properties of Frobenius norm)

We have the following properties.

  1. Let A , B R m × m be symmetric and skew-symmetric, respectively, i.e. A = A T , B = - B T . Then

    (4.2) A + B F 2 = A F 2 + B F 2 .

  2. The Frobenius norm is unitary invariant, i.e.

    (4.3) A F = U A V F

    for A R m × m and U , V R m × m unitary.

  3. If A SPD with (positive) eigenvalues ( λ i ) i = 1 m , then

    (4.4) Log ( A ) F = ( i = 1 m log 2 ( λ i ) ) 1 / 2 .

Proof

The proof of the first item is straightforward by using the definition of F in (2.4). The second item follows directly by considering the trace representation of the Frobenius norm [51]:

U A V F 2 = trace ( ( U A V ) T U A V ) = trace ( V T A T A V ) = trace ( A V V T A T ) = A F 2 .

The third item is a direct consequence of Remark 4.1(i), Definition 4.2(ii) and (4.3). ∎

The last lemma of this subsection deals with the set SPD [ ε ¯ , ε ¯ ] spec , the set of symmetric, positive definite matrices with bounded spectrum in the interval [ ε ¯ , ε ¯ ] , defined in (2.2). We need this result later in the numerical implementation for defining a suitable projection.

Given an arbitrary matrix A R m × m , there always exists a unique matrix M SPD [ ε ¯ , ε ¯ ] spec which is closest in the Frobenius norm, i.e.

M = arg min X SPD [ ε ¯ , ε ¯ ] spec A - X F 2 .

The minimizing matrix 𝑀 can be computed explicitly as stated in the following lemma. The proof is done in a way similar to [25, Theorem 2.1] and is included here for completeness.

Lemma 4.5

Let A R m × m . Define B := 1 2 ( A + A T ) and C := 1 2 ( A - A T ) as the symmetric and skew-symmetric parts of 𝐴, respectively. Let ( λ i ) i = 1 m be the eigenvalues of 𝐵 which can be decomposed into B = Z Λ Z T , where 𝑍 is a unitary matrix, i.e. Z Z T = Z T Z = 1 m , and Λ = diag ( λ 1 , , λ m ) . Then the unique minimizer of

(4.5) min X SPD [ ε ¯ , ε ¯ ] spec A - X F 2 ,

where SPD [ ε ¯ , ε ¯ ] spec is defined in (2.2), is Z Y Z T with Y = diag ( d 1 , , d m ) , where

d i := { λ i if λ i [ ε ¯ , ε ¯ ] , ε ¯ if λ i > ε ¯ , ε ¯ if λ i < ε ¯ .

Proof

By definition of 𝐵 and 𝐶, we can write A = B + C and thus

A - X F 2 = B + C - X F 2 = B - X F 2 + C F 2 ,

where we used (4.2) in the second equality. The problem in (4.5) thus reduces to finding

arg min X SPD [ ε ¯ , ε ¯ ] spec B - X F 2 .

The matrix 𝐵 is symmetric, and thus we can write B = Z Λ Z T , where Z R m × m is a unitary matrix whose columns are the eigenvectors of 𝐵 and Λ R m × m is a diagonal matrix whose entries are the eigenvalues of 𝐵, i.e. Λ = diag ( λ 1 , , λ m ) . Let Y = Z T X Z be similar to 𝑋 so that spec ( Y ) [ ε ¯ , ε ¯ ] . Then we obtain, by using (4.3),

B - X F 2 = Λ - Y F 2 = { i , j : i j } y i j 2 + i = 1 m ( λ i - y i i ) 2 = { i , j : i j } y i j 2 + { i : λ i [ ε ¯ , ε ¯ ] } ( λ i - y i i ) 2 + { i : λ i > ε ¯ } ( λ i - y i i ) 2 + { i : λ i < ε ¯ } ( λ i - y i i ) 2 { i : λ i > ε ¯ } ( λ i - ε ¯ ) 2 + { i : λ i < ε ¯ } ( λ i - ε ¯ ) 2 .

Thus the lower bound is uniquely attained for Y := diag ( d i ) with

d i := { λ i if λ i [ ε ¯ , ε ¯ ] , ε ¯ if λ i > ε ¯ , ε ¯ if λ i < ε ¯ .

4.2 Existence

After giving the needed definitions from matrix calculus, the goal of this subsection is now to apply Theorem 3.3 to the regularization functional defined in (1.2) with the set K = SPD z Log defined in (4.1) and associated log-Euclidean metric defined below in (4.6). Therefore, we need to prove the equivalence of the log-Euclidean and Euclidean metric to guarantee in particular that Assumption 2.1(v) is fulfilled. Then Assumption 3.1 holds true as stated in Remark 3.2, and therefore Theorem 3.3 is applicable.

We start by defining and reviewing some properties of the log-Euclidean metric.

Definition 4.6

Definition 4.6 (Log-Euclidean metric)

Let A , B SPD . The log-Euclidean metric is defined as

(4.6) d SPD ( A , B ) := d ( A , B ) := Log ( A ) - Log ( B ) F , A , B SPD .

Lemma 4.7

The log-Euclidean metric satisfies the metric axioms on SPD .

Proof

This follows directly because F is a norm and Log restricted to SPD is a diffeomorphism. ∎

The reasons for choosing this measure of distance is stated in the following remark.

Remark 4.8

The log-Euclidean metric arises when considering SPD not just as convex cone in the vector space of matrices but as a Riemannian manifold. Thus it can be endowed with a Riemannian metric defined by an inner product on the tangent space, see for example [18, 38, 3]. Two widely used geodesic distances are the affine-invariant metric

(4.7) d AI ( A , B ) = Log ( A - 1 / 2 B A - 1 / 2 ) F , A , B SPD ,

and the log-Euclidean metric as stated above. These measures of dissimilarity are more adequate in DTMRI as pointed out in [3] because zero or negative eigenvalues induce an infinite distance.

The affine-invariant distance measure is computationally much more demanding, which is a major drawback. This is not the case for the log-Euclidean distance, which leads to Euclidean distance computations in the matrix logarithmic domain.

The following statement can be found in [33, Sections 2.4 & 2.5].

Lemma 4.9

Let d = d SPD denote the log-Euclidean metric (as defined in (4.6)), d AI the affine-invariant metric (as defined in (4.7)) and d R m × m the standard Euclidean distance. Then ( SPD , d ) as well as ( SPD , d AI ) form a complete metric space. This is not the case for ( SPD , d R m × m ) .

Different metrics induce different properties on the corresponding regularizer. We compare Φ [ d SPD ] l , Φ [ d R m × m ] l and additionally a Sobolev semi-norm regularizer in the following.

Remark 4.10

Remark 4.10 (Invariances)

For the log-Euclidean metric d = d SPD as defined in (4.6), the following holds true.

  • Scale invariance: Let c > 0 and A , B SPD , and denote by 1 m the identity matrix in R m × m . Then

    d ( c A , c B ) = d ( c 1 m A , c 1 m B ) = Log ( c 1 m ) + log ( A ) - Log ( c 1 m ) - Log ( B ) F = d ( A , B ) .

  • Invariance under inversion: Let A , B SPD . Because Log ( A - 1 ) = - Log ( A ) , we directly get that

    d ( A - 1 , B - 1 ) = d ( A , B ) .

  • Unitary invariance: Let A , B SPD and 𝑈 unitary. Because of the unitary invariance of the Frobenius norm, then

    d ( U A U T , U B U T ) = d ( A , B ) .

These properties transfer to our regularizer Φ [ d ] l over W s , p ( Ω ; SPD ) . Clearly, when considering Φ [ d R m × m ] l , where d R m × m ( A , B ) = A - B F is the standard Euclidean distance, the first two properties do not hold true in contrast to the unitary invariance which is also valid.

Although we only work with fractional derivatives of order s ( 0 , 1 ) , we consider for comparison purposes the regularization functional (see also (6.1))

w W 1 , p ( Ω , R m × m ) Θ ( w ) := Ω w ( x ) F p d x .

None of the invariances above, i.e. scale invariance, invariance under inversion and unitary invariance, is valid for Θ.

Instead,

Θ ( w + C ) = Θ ( w ) , Θ ( w ) = Θ ( - w )

for some constant matrix C R m × m , i.e. it is translation and reflection invariant. This, in turn, does not hold (or is not even well-defined) for our regularizer Φ [ d ] l with the log-Euclidean metric but as well when considering the standard Euclidean distance, i.e. it does hold for Φ [ d R m × m ] l . A comparison is shown in Table 1.

Table 1

Comparison of invariance properties of our regularizer Φ [ d ] l , Φ [ d R m × m ] l and the regularization term Θ.

Φ [ d ] l Φ [ d R m × m ] l Θ
scale invariant
inversion invariant
unitary invariant
translation invariant
reflection invariant

In order to show that Theorem 3.3 is applicable for F [ d ] α , w δ defined in (1.2) with K = SPD z Log and associated log-Euclidean metric d = d SPD defined in (4.6), we have to show that Assumption 3.1 and therefore Assumption 2.1, in particular the equivalence stated in (v), is valid. In order to prove that, we need the following corollary.

Corollary 4.11

Let A SPD z Log (defined in (2.3)) with eigenvalues ( λ i ) i = 1 m . Then, for each i = 1 , , m ,

λ i [ e - z , e z ] ,

i.e. SPD z Log SPD [ e - z , e z ] spec (for the definition of the latter set, see (2.2)).

Proof

If A SPD z Log , it holds that Log ( A ) F z . Using (4.4), this is equivalent to i = 1 m log 2 ( λ i ) z 2 , so the claim follows. ∎

Note that the reverse embedding in the previous lemma does not hold true. If A SPD [ e - z , e z ] spec such that, for each eigenvalue λ i , i = 1 , , m , we have that λ i { e - z , e z } , then A SPD z m Log SPD z Log .

Now we can prove that the Euclidean and the log-Euclidean metric are equivalent on SPD z Log . In particular, we show that Log is bi-Lipschitz on SPD z Log and calculate the constant explicitly. Without explicit computation, this would follow from the fact that Log is a diffeomorphism on symmetric, positive definite matrices and that SPD z Log is a compact subset.

Lemma 4.12

Let A , B SPD z Log be defined in (2.3). Then

(4.8) 1 e z A - B F 2 Log ( A ) - Log ( B ) F 2 1 e - z A - B F 2 .

Proof

Since 𝐴 and 𝐵 are symmetric and positive definite, they can be factorized using their eigendecomposition, see Remark 4.1(i). Hence we can write

A = U Λ A U T , B = V Λ B V T ,

where U , V R m × m are unitary matrices and Λ A , Λ B are diagonal matrices whose entries are the corresponding positive eigenvalues ( λ 1 , , λ m ) of 𝐴 and ( μ 1 , , μ m ) of 𝐵, respectively. By Corollary 4.11, it holds that λ i , μ i [ e - z , e z ] for all i = 1 , , m . We consider two cases.

Case 1. We assume that all eigenvalues of 𝐴 and 𝐵 are equal, i.e. they have the same one-dimensional spectrum spec ( A ) = spec ( B ) = { λ } , meaning that Λ := λ 1 m := Λ A = Λ B . This in turn gives that

A - B F 2 = U Λ U T - V Λ V T F 2 = V T U Λ - Λ V T U F 2 = V T U λ - λ V T U F 2 = 0 , Log ( A ) - Log ( B ) F 2 = V T U log ( Λ ) - log ( Λ ) V T U F 2 = V T U log ( λ ) - log ( λ ) V T U F 2 = 0 ,

using the unitary invariance of the Frobenius norm as stated in (4.3) and the properties of the matrix logarithm in Definition 4.2(ii) in the second equation. Thus (4.8) is trivially fulfilled.

Case 2. We now assume that there exist at least two different eigenvalues λ i μ j , i , j { 1 , , m } of 𝐴 and 𝐵.

We show the lower inequality 1 e z A - B F 2 Log ( A ) - Log ( B ) F 2 in (4.8). The upper inequality can be done analogously.

By (4.3) and the properties of the matrix logarithm in Definition 4.2(ii), it follows that

(4.9) Log ( A ) - Log ( B ) F 2 = U Log ( Λ A ) U T - V Log ( Λ B ) V T F 2 = V T U Log ( Λ A ) - Log ( Λ B ) V T U F 2 = C ( diag ( log ( λ 1 ) , , log ( λ m ) ) ) - diag ( log ( μ 1 ) , , log ( μ m ) ) C F 2 ,

where C := V T U . Using the definition of the Frobenius norm in (2.4), we obtain further that

(4.10) C ( diag ( log ( λ 1 ) , , log ( λ m ) ) ) - diag ( log ( μ 1 ) , , log ( μ m ) ) C F 2 = i , j = 1 m | c i j ( log ( λ j ) - log ( μ i ) ) | 2 .

Indices ( i , j ) { 1 , , m } for which λ j = μ i do not contribute to the sum in (4.10) (and do not change the following calculation), so we define I := { ( i , j ) { 1 , , m } : λ j μ i } as the set of such indices ( i , j ) { 1 , , m } for which we have λ j μ i .

From the mean value theorem, it follows that, for every ( i , j ) I , there exists some

ξ i j { ( λ j , μ i ) [ e - z , e z ] if λ j < μ i , ( μ i , λ j ) [ e - z , e z ] if μ i < λ j ,

such that

(4.11) ( i , j ) I | c i j ( log ( λ j ) - log ( μ i ) ) | 2 = ( i , j ) I | c i j 1 ξ i j ( λ j - μ i ) | 2 1 e z ( i , j ) I | c i j ( λ j - μ i ) | 2 .

Further, we can write

(4.12) 1 e z ( i , j ) I | c i j ( λ j - μ i ) | 2 = C ( diag ( λ 1 , , λ m ) ) - diag ( μ 1 , , μ m ) C F 2 = 1 e z C Λ A - Λ B C F 2 .

Combining (4.9), (4.10), (4.11), (4.12), the definition of C = V T U and (4.3), we obtain that

Log ( A ) - Log ( B ) F 2 1 e z C Λ A - Λ B C F 2 = 1 e z U Λ A U T - V Λ B V T F 2 = 1 e z A - B F 2 ,

which finishes the proof. ∎

The previous lemma, Lemma 4.12, proves that Assumption 2.1(v) is valid. This together with Remark 3.2 proves the following theorem.

Theorem 4.13

Let K = SPD z Log and d = d SPD as in (4.6). Then the functional F [ d ] α , w δ as defined in (1.2) over W s , p ( Ω ; SPD z Log ) satisfies the assertions of Theorem 3.3. In particular, it attains a minimizer and fulfills a stability and convergence result.

4.3 Uniqueness

So far we showed that the functional F [ d ] α , w δ as defined in (1.2) over W s , p ( Ω ; SPD z Log ) using the log-Euclidean metric d = d SPD as in (4.6) attains a minimizer. In this subsection, we prove that the minimum is unique. To this end, we consider the symmetric, positive definite matrices from a differential geometric point of view.

The following lemma can be found in [3] and also [39].

Lemma 4.14

The space ( SPD , d ) , where d = d SPD denotes the log-Euclidean metric as defined in (4.6), is a complete, connected Riemannian manifold with zero sectional curvature.

In other words, ( SPD , d ) is a flat Hadamard manifold and therefore in particular a Hadamard space. The last property guarantees that the metric 𝑑 is geodesically convex [48, Corollary 2.5], i.e. let γ , η : [ 0 , 1 ] SPD be two geodesics; then

(4.13) d ( γ t , η t ) t d ( γ 0 , η 0 ) + ( 1 - t ) d ( γ 1 , η 1 ) .

Moreover, d p is strictly convex in one argument for p > 1 ([48, Proposition 2.3] & [4, Example 2.2.4]), i.e. for M SPD fixed and γ 0 γ 1 ,

(4.14) d p ( γ t , M ) < t d p ( γ 0 , M ) + ( 1 - t ) d p ( γ 1 , M ) .

The following result states that connecting geodesics between two points in SPD z Log stay in this set.

Lemma 4.15

Let Assumption 3.1 hold. Let K = SPD z Log , and let d = d SPD be the log-Euclidean metric as defined in (4.6). Let w * , w W s , p ( Ω ; SPD z Log ) W s , p ( Ω ; SPD ) . For γ : Ω × [ 0 , 1 ] SPD , define

γ x := γ ( x , ) : [ 0 , 1 ] W s , p ( Ω ; SPD ) ,

as a connecting geodesic between γ x ( 0 ) = w * ( x ) , and γ x ( 1 ) = w ( x ) and

γ t := γ ( , t ) : Ω W s , p ( Ω ; SPD ) ,

as the evaluation of the geodesic between w * ( x ) and w ( x ) at time 𝑡 for x Ω . Then γ t W s , p ( Ω ; SPD z Log ) .

Proof

We split the proof into two parts. First we show that γ t maps indeed into SPD z Log . Afterwards, we prove that it actually lies in W s , p ( Ω ; SPD z Log ) .

The geodesic γ t connects γ 0 ( x ) = w * ( x ) and γ 1 ( x ) = w ( x ) for x Ω . Therefore ([47, Chapter 3.5] and [3]), it can be written as

γ t ( x ) = Exp ( t Log ( w * ( x ) ) + ( 1 - t ) Log ( w ( x ) ) )

which is equivalent to

Log ( γ t ( x ) ) = t Log ( w * ( x ) ) + ( 1 - t ) Log ( w ( x ) ) .

We denote by 1 m the identity matrix of size m × m and note that Log ( γ t ( x ) ) F = d ( γ 0 ( x ) , 1 m ) , where the Log-Euclidean metric 𝑑 is as defined in (4.6). Because of the geodesic convexity, see (4.13), we obtain that

Log ( γ t ( x ) ) F = d ( γ t ( x ) , 1 m ) t d ( γ 0 ( x ) , 1 m ) + ( 1 - t ) d ( γ 1 ( x ) , 1 m ) t z + ( 1 - t ) z = z

because w * , w W s , p ( Ω ; SPD z Log ) , i.e. Log ( w * ) F z , Log ( w ) F z . This shows that γ t maps into SPD z Log .

Next we need to prove that actually γ t W s , p ( Ω ; SPD z Log ) , i.e. that

γ t W s , p ( Ω ; R m × m ) p = Ω γ t ( x ) F p d x + Ω × Ω γ t ( x ) - γ t ( y ) F p | x - y | n + p s d ( x , y ) = Ω γ t ( x ) F p d x + Φ [ d R m × m ] 0 ( γ t ) < .

We obtain by Jensen’s inequality that

γ t W s , p ( Ω ; R m × m ) p 2 p - 1 ( Ω γ t ( x ) - 1 m F p d x + Ω 1 m F p d x ) + Φ [ d R m × m ] 0 ( γ t ) .

Using (4.8), it follows that

2 p - 1 ( Ω γ t ( x ) - 1 m F p d x + Ω 1 m F p d x ) + Φ [ d R m × m ] 0 ( γ t ) 2 p - 1 ( e z ) p / 2 ( Ω d p ( γ t ( x ) , 1 m ) d x + Φ [ d ] 0 ( γ t ) ) + C ,

where C := 2 p - 1 | Ω | . By using the geodesic convexity stated in (4.13) and (4.14) and again the equivalence of the Euclidean and the log-Euclidean metric (see Lemma 4.12), we get that

2 p - 1 ( e z ) p / 2 ( Ω d p ( γ t ( x ) , 1 m ) d x + Φ [ d ] 0 ( γ t ) ) + C 2 p - 1 ( e z ) p / 2 ( t Ω d p ( γ 0 ( x ) , 1 m ) d x + ( 1 - t ) Ω d p ( γ 1 ( x ) , 1 m ) d x + t Φ [ d ] 0 ( γ 0 ) + ( 1 - t ) Φ [ d ] 0 ( γ 1 ) ) + C 2 p - 1 e p z ( t Ω γ 0 ( x ) - 1 m F p d x + ( 1 - t ) Ω γ 1 ( x ) - 1 m F p d x + t Φ [ d R m × m ] 0 ( γ 0 ) + ( 1 - t ) Φ [ d R m × m ] 0 ( γ 1 ) ) + C .

The last expression is finite because of the assumption that w * , w W s , p ( Ω ; SPD z Log ) . ∎

Now we can state the uniqueness result.

Theorem 4.16

Let Assumption 3.1 hold. Let K = SPD z Log , and let d = d SPD be the log-Euclidean metric as defined in (4.6). Then the functional F [ d ] α , w δ as defined in (1.2) on W s , p ( Ω ; SPD z Log ) attains a unique minimizer.

Proof

Existence of a minimizer is guaranteed by Theorem 4.13. Now, let us assume that there exist two minimizers w * w W s , p ( Ω ; SPD z Log ) of the functional F [ d ] α , w δ .

Analogously as in Lemma 4.15 for a geodesic path γ : Ω × [ 0 , 1 ] SPD connecting w * and w , we denote γ t = γ ( , t ) for t [ 0 , 1 ] . Thus, in particular, w * ( x ) = γ 0 ( x ) and w ( x ) = γ 1 ( x ) for x Ω . Especially, γ t W s , p ( Ω ; SPD z Log ) (see Lemma 4.15).

Because w δ is fixed, 𝑑 is strictly convex in one argument by (4.14) and convex in both arguments by (4.13), it follows that

(4.15) F [ d ] α , w δ ( γ t ) = Ω χ Ω D ( x ) d p ( γ t ( x ) , w δ ( x ) ) d x + α Ω × Ω d p ( γ t ( x ) , γ t ( y ) ) | x - y | n + p s ρ l ( x - y ) d ( x , y ) < t F [ d ] α , w δ ( γ 0 ) + ( 1 - t ) F [ d ] α , w δ ( γ 1 ) .

Because w * and w are both minimizers, we have that

F [ d ] α , w δ ( w * ) = F [ d ] α , w δ ( γ 0 ) = F [ d ] α , w δ ( γ 1 ) = F [ d ] α , w δ ( w ) .

In particular, for t = 1 / 2 , we obtain by the above equality and by (4.15) that

F [ d ] α , w δ ( γ 1 / 2 ) < 1 2 F [ d ] α , w δ ( γ 0 ) + 1 2 F [ d ] α , w δ ( γ 1 ) = F [ d ] α , w δ ( γ 0 ) = min w W s , p ( Ω ; SPD z Log ) F [ d ] α , w δ ( w ) ,

which is a contradiction to the minimizing property of w * ( x ) = γ 0 ( x ) and w ( x ) = γ 1 ( x ) for x Ω . Hence γ 0 and γ 1 must be equal forcing equality in (4.15) and thus giving that the minimum is unique. ∎

Existence and uniqueness in the case s p > n

If s p > n , then existence and uniqueness of the minimizer of the functional F [ d ] α , w δ even holds on the larger set W s , p ( Ω ; SPD ) rather than on W s , p ( Ω ; SPD z Log ) , where SPD is associated with the log-Euclidean distance d = d SPD as defined in (4.6). Existence in Theorem 4.13 and uniqueness in Theorem 4.16 (with K = SPD z Log ) are based on the theory provided in [14] (see Theorem 3.3), where it is a necessary assumption that the set 𝐾 is closed which is not the case for the set SPD .

Nevertheless, it is possible to get existence and uniqueness on this set because, for every minimizing sequence w k W s , p ( Ω ; SPD ) , k N , we automatically get that w k W s , p ( Ω ; SPD z Log ) so that it takes values on the closed subset SPD z Log . Then existence of a unique minimizer on W s , p ( Ω ; SPD ) follows by the proofs already given, see [14, Theorem 3.6] and Theorem 4.16.

We now sketch the proof of the assertion. Throughout this sketch, 𝐶 denotes a finite generic constant which, however, can be different from line to line.

Sketch of assertion

Denote by d = d SPD the log-Euclidean metric (as defined in (4.6)). Let us take a minimizing sequence w k W s , p ( Ω ; SPD ) , k N , of F [ d ] α , w δ so that we can assume that F [ d ] α , w δ ( w k ) C for all w k , k k 0 N .

Computing the log-Euclidean metric leads to evaluations of the Euclidean metric in the matrix logarithmic domain, cf. (4.6), meaning that

d ( A , B ) = Log ( A ) - Log ( B ) F = d R m × m ( Log ( A ) , Log ( B ) ) , A , B SPD .

Because of this and the fact that w δ L p ( Ω ; SPD ) , we get that

C Log ( w k ) L p ( Ω ; SYM ) p + α Φ [ d R m × m ] 1 ( Log ( w k ) ) .

Because of [32, Lemma 2.7], we can thus bound the W s , p -norm of Log ( w k ) ,

(4.16) C Log ( w k ) W s , p ( Ω ; SYM ) p .

If s p > n , the space W s , p ( Ω ; R m × m ) is embedded into Hölder spaces C 0 , α ( Ω ; R m × m ) with α := ( s p - n ) / p guaranteed by [17, Theorem 8.2]. Because of (4.16), this gives us that

C Log ( w k ) C 0 , α ( Ω ; SYM ) p ,

yielding in particular that Log ( w k ) < C := z .

By the definition of SPD z Log in (2.3), we thus obtain that w k W s , p ( Ω ; SPD z Log ) for all k k 0 . Hence every minimizing sequence w k W s , p ( Ω ; SPD ) , k N , of F [ d ] α , w δ is automatically a minimizing sequence in W s , p ( Ω ; SPD z Log ) .

5 Numerics

In this section, we go into more detail on the minimization of the regularization functional F [ d ] α , w δ defined in (1.2) with the log-Euclidean metric d = d SPD as defined in (4.6) (see [3]) over the set K = SPD z Log , the set of symmetric, positive definite m × m matrices with bounded logarithm, as defined in (4.1) for denoising and inpainting of DTMRI images.

To optimize F [ d ] α , w δ , we use a projected gradient descent algorithm. The implementation is done in Matlab. The gradient step is performed by using Matlab’s built-in function fminunc, where the gradient is approximated with a finite difference scheme (central differences in the interior and one-sided differences at the boundary). Therefore, after each step, we project the data which are elements of the larger space SYM ( 3 ) back onto K = SPD z Log ( 3 ) by applying the following projection. We remark that, by projection, we here refer to an idempotent mapping.

5.1 Projections

Definition 5.1

Definition 5.1 (Projection operators)

We define the following projections.

  • Projection of SYM onto SPD [ ε , ) spec : Let M SYM be a symmetric matrix with eigendecomposition M = V Λ V T with Λ = diag ( λ 1 , , λ m ) . Then the projection of 𝑀 onto the set SPD [ ε , ) spec is given by

    (5.1) P 1 : SYM SPD [ ε , ) spec , M V Σ V T ,

    where Σ = diag ( μ 1 , , μ m ) with

    μ i := { λ i if λ i ε , ε if λ i < ε .

  • Projection of SPD [ ε , ) spec onto SPD z Log : Let M SPD [ ε , ) spec with eigenvalues ( λ i ) i = 1 m and eigendecomposition M = V Λ V T . Define

    C Frob := Log ( M ) F 2 = i = 1 m log 2 ( λ i )

    as the squared Frobenius norm of Log ( M ) . Then the pr