Matrix-Sequences of Geometric Means in the Case of Hidden (Asymptotic) Structures

1. Introduction

The study of geometric means of matrices has become increasingly important in various fields, including numerical linear algebra, differential equations, and signal processing. Positive definite matrices arise naturally in many applications, such as in the discretization of differential equations via methods like finite differences or finite elements. These methods produce matrix-sequences whose sizes increase as the mesh size becomes finer. An essential concept related to matrix-sequences is the asymptotic eigenvalue distribution in the Weyl sense and this has been a classical study; see e.g. [17,23,40]. Since the publication of the seminal papers by Tyrtyshnikov in 1996 [41] and Tilli in 1998 [38,39], there has been growing interest in this area, which eventually contributed to the development of the theory of Generalized Locally Toeplitz (GLT) sequences [35,36]. The increasing attention to this topic is not purely theoretical, as the asymptotic eigenvalue and singular value distributions have significant practical applications, particularly in the analysis of large-scale matrix computations especially in the context of the numerical approximations of systems of (fractional) partial differential equations ((F)PDEs); see e.g. the books and review papers [7,8,20,21,22] and references therein. In fact, it is worth noticing that virtually any meaningful approximation technique of a (F)PDE leads to a GLT matrix-sequence, including finite differences, finite elements of any order, discontinuous Galerkin techniques, finite volumes, isogeometric analysis etc. More in detail, if we fix $d,r$ positive integer numbers, then the set of d-level r-block GLT matrix-sequences forms a maximal *-algebra of matrix-sequences, isometrically equivalent to d-variate $r\times r$ matrix-valued mesurable functions defined on ${[0,1]}^d\times {[-\pi ,\pi ]}^d$. Furthermore, a d-level r-block GLT sequence ${\left\{A_n\right\}}_n$ is uniquely associated with a $r\times r$ matrix-valued Lebesgue-measurable function $\kappa$, known as the GLT symbol, which is defined over the domain $D={[0,1]}^d\times {[-\pi ,\pi ]}^d$. Notice that the set ${[0,1]}^d$ can be replaced by any bounded Peano-Jordan measurable subset of ${\mathbb{R}}^d$ as occurring with the notion of reduced GLT *-algebras, see [35][pp. 398-399, formula (59)] for the first occurrence with applications of approximated PDEs on general non-Cartesian domains in d dimensions, [36][Section 3.1.4] for the first formal proposal, and [5] for an exhaustive treatment, containing both the *-algebra theoretical results and a number of applications. This symbol provides a powerful tool for analyzing the singular value and eigenvalue distributions when the matrices $A_n$ are Hermitian and part of a matrix-sequence of increasing size. The notation ${\left\{A_n\right\}}_n{\sim }_{\mathrm{GLT}}\kappa$ indicates that ${\left\{A_n\right\}}_n$ is a GLT sequence with symbol $\kappa$. Notably, the symbol of a GLT sequence is unique in the sense that if ${\left\{A_n\right\}}_n{\sim }_{\mathrm{GLT}}\kappa$ and ${\left\{A_n\right\}}_n{\sim }_{\mathrm{GLT}}\xi$, then $\kappa =\xi$ almost everywhere in ${[0,1]}^d\times {[-\pi ,\pi ]}^d$[7,8,20,21]. Furthermore, by the *-algebra structure, ${\left\{A_n\right\}}_n{\sim }_{\mathrm{GLT}}\kappa$ and ${\left\{B_n\right\}}_n{\sim }_{\mathrm{GLT}}\kappa$ implies that ${\{A_n-B_n\}}_n{\sim }_{\mathrm{GLT}}0$, i.e., the matrix-sequence ${\{A_n-B_n\}}_n$ is zero-distributed and the latter is very important for building explicit matrix-sequences approximating a GLT matrix-sequence and whose inversion is computationally cheap, in the context of preconditioning of large linear systems.

In certain physical applications, it is often necessary to represent the results of multiple experiments through a single average matrix G, where the data is represented by a set of positive definite matrices $A_1,A_2,...,A_k$. The arithmetic mean ${\displaystyle \frac{1}{k}}{\sum }_{i=1}^k{A}_i$ is not appropriate in these cases because it does not fulfill the requirement that the inverse of the mean should coincide with the average of the inverses $G^{-1}={\displaystyle \frac{1}{k}}{\sum }_{i=1}^k{A}_i^{-1}$. This property, which is crucial in certain physical models is satisfied by the geometric mean. For positive real numbers $a_1,a_2,...,a_k$, the geometric mean is defined as

$$g={\left(\prod _{i=1}^k{a}_i\right)}^{1/k},$$

a concept that is extended to the case of matrices [10,26] in a nontrivial way, where the difficulty is of course the lack of commutativity. A well-known definition that satisfies desirable properties such as congruence invariance, permutation invariance, and consistency with the scalar geometric mean was proposed by ALM [1]. They defined the geometric mean of two Hermitian positive definite (HPD) matrices as

$$G(A,B)=A^{1/2}{\left(A^{-1/2}B{A}^{-1/2}\right)}^{1/2}{A}^{1/2}.$$

(1)

We recall the definition of functions applied to diagonalizable matrices, which we frequently utilize throughout our work. Suppose $A\in {\mathbb{C}}^{n\times n}$ is diagonalizable, meaning $A=MD{M}^{-1}$, where $D=\left(d_{ij}\right)$ is a diagonal matrix, and f is a given function. In this case, $f\left(A\right)$ is defined as $Mf\left(D\right)M^{-1}$, with $f\left(D\right)$ being a diagonal matrix whose diagonal entries are $f\left(d_{ii}\right)$, for $i=1,\cdots ,n$. In the case where f is a multi-valued function, the same branch of f must be chosen for any repeated eigenvalue.

If $k>2$ then the ALM-mean is obtained through a recursive iteration process where at each step the geometric means of k matrices are computed by reducing them to $k-1$ matrices. However, a significant limitation of this method is its linear convergence leading to a high computational cost due to the large number of iterations required at each recursive step. As a result, the computation of the ALM-mean using this approach becomes quite expensive. However, despite the elegance of the ALM geometric mean for two matrices, it becomes computationally infeasible to extend this formula to more than two matrices [14].

To overcome these limitations, the Karcher mean [13], was introduced as a generalization of the geometric mean for more than two matrices. The Karcher mean of HPD matrices $A_1,A_2,\cdots ,A_k$ is defined as the unique positive definite solution to the matrix equation:

$$\sum _{i=1}^{k}log\left(A_i^{-1}X\right)=0,$$

(2)

as established by Moakher [26, Proposition 3.4]. This equation can be equivalently expressed in other forms, such as

$$\sum _{i=1}^{k}log\left(X{A}_i^{-1}\right)=0,\mathrm{or}\sum _{i=1}^{k}log\left(X^{1/2}{A}_i^{-1}{X}^{1/2}\right)=0,$$

(3)

by utilizing the formula $M^{-1}log\left(K\right)M=log\left(M^{-1}KM\right)$, which holds for any invertible matrix M and any matrix K with real positive eigenvalues. This formulation arises from Riemannian geometry, where the space of positive definite matrices forms a Riemannian manifold with non-positive curvature. The Karcher mean represents the center of mass (or barycenter) on this manifold [10]. In this manifold, the distance between two positive definite matrices A and B is defined as $\delta (A,B)=\parallel log\left(A^{-1/2}B{A}^{-1/2}\right){\parallel }_F,$ where ${\parallel ·\parallel }_F$ denotes the Frobenius norm.

Several numerical methods have been proposed for solving the Karcher mean equation. Initially, fixed-point iteration methods were used, but these methods suffered from slow convergence, especially in cases where the matrices involved were poorly conditioned [26]. Later, methods based on gradient descent in the Riemannian manifold were introduced. A common iteration scheme for approximating the Karcher mean is

$$X_{v+1}=X_{v}exp\left(-{\theta }_v\sum _{i=1}^{k}log\left(A_i^{-1}{X}_v\right)\right),X_0\in {\mathcal{P}}_n,v\ge 0,$$

(4)

where $X_v$ is the approximation at the v-th step, and the exponential and logarithmic functions are matrix operations. Although this method improves convergence, it can exhibit slow linear convergence in certain cases. The iteration with $X_0=A_1$ or $X_0=I$ and ${\theta }_v=1/k$, as considered in [25,26], can fail to converge for some matrices $A_1,\cdots ,A_k$. Furthermore, similar iterations have been proposed in [3,28], but without specific recommendations on choosing the initial value or step length. While an optimal ${\theta }_v$ could theoretically be determined using a line search strategy, this approach is often computationally expensive. Heuristic strategies for selecting the step size, as discussed in [18], may result in slow convergence in many cases.

To further enhance the convergence rate, the Richardson iteration method is employed. Indeed, the considered method improves the convergence by using a parameter $\theta$, which controls the step size in each iteration [13]. More precisely, given $X_0\in {\mathcal{P}}_n$, the Richardson iteration is given by

$$X_{v+1}=X_v-\theta X_v\sum _{i=1}^{k}log\left(A_i^{-1}{X}_v\right)\equiv T\left(X_v\right),v\ge 0.$$

(5)

Any solution of Eq. (3) is a fixed point of the map T in (5). The iterative formula can also be rewritten as

$$X_{\nu +1}=X_{\nu }-\theta X_{\nu }^{1/2}\left(\sum _{i=1}^{k}log\left(X_{\nu }^{1/2}{A}_i^{-1}{X}_{\nu }^{1/2}\right)\right)X_{\nu }^{1/2},v\ge 0,$$

(6)

provided that all the iterates $X_{\nu }$ remain positive definite. Equation (6) further demonstrates that if $X_{\nu }$ is Hermitian, then $X_{\nu +1}$ is also a Hermitian matrix.

The parameter $\theta$ plays a crucial role in controlling the step size of each iteration and choosing an optimal value of $\theta$ can significantly influence the convergence behavior of the iteration. If $\theta$ is small enough, the iteration is guaranteed to produce positive definite matrices and converge towards the solution. In particular, when the matrices $A_1,\cdots ,A_k$ commute, setting $\theta ={\displaystyle \frac{1}{k}}$ ensures at least quadratic convergence. More generally, the optimal value of $\theta$ can be determined based on the condition numbers of matrices $M_i=G^{1/2}{A}_i^{-1}{G}^{1/2},i=1,\cdots ,k$, where G is the desired solution. The closer the eigenvalues of $M_i$ are to 1, the faster the convergence is. This analysis guarantees that if the initial guess $X_0$ is close enough to the solution and is positive definite, the sequence $\left\{X_v\right\}$ generated by the iteration remains well-defined and converges to the desired solution. However, if the initial iterate is not positive definite, adjusting the value of $\theta$ or modifying the iteration scheme may be necessary to ensure that all iterates remain positive.

Numerical experiments show that selecting appropriate initial guesses, such as the arithmetic mean $X_0={\displaystyle \frac{1}{k}}{\sum }_{i=1}^k{A}_i$ or the identity matrix $X_0=I$, can significantly affect the convergence rate. In particular, the cheap mean introduced in [12] provides a practical initial approximation that leads to faster convergence in many cases.

In our study, we consider the analysis of the Karcher mean to matrix-sequences, particularly those arising from the discretization of differential equations, which often form GLT sequences. Numerical results demonstrate that the geometric mean of not only two but more than two GLT matrix-sequences is itself a GLT matrix-sequence, with the symbol of the new sequence given by the geometric mean of the original symbols: the latter is formally proven in the case of two GLT matrix-sequences. Regarding the examples, we consider either scalar unilevel and multilevel GLT matrix-sequences or block GLT asymptotic structures, with special attention to cases stemming from the approximation by local methods of differential operators. By analyzing the spectral distribution of the geometric mean of these matrix-sequences, we provide new insights into the asymptotic behavior of large-scale matrix computations and their potential applications in numerical analysis.

This paper is structured as follows. In Section 2, we introduce notations, terminology, and preliminary results concerning Toeplitz and GLT structures, which are essential for the mathematical formulation of the problem and its technical solution. In Section 3, we present the geometric mean of two matrices and the Karcher mean for more than two matrices, followed by a discussion of the iterative methods employed for their computation: the section contains the GLT theoretical results. Section 4 contains numerical experiments that illustrate the (asymptotic) spectral behavior of the geometric mean for GLT matrix-sequences in both 1D and 2D settings and in both scalar and block cases. Finally, in Section 5 we draw conclusions and we punt in evidence few open problems.

2. Preliminaries

In this introduction, we provide the necessary tools for performing the spectral analysis of the matrices involved, based on the theory of uni and multilevel block GLT matrix-sequences.

2.1. Matrices and Matrix-Sequences

Given a square matrix $A\in {\mathbb{C}}^{m\times m}$, we denote by $A^{*}$ its conjugate transpose and by $A^{†}$ the Moore–Penrose pseudoinverse of A. Recall that $A^{†}=A^{-1}$ whenever A is invertible. Regarding matrix norms $\parallel ·\parallel$ refers to the spectral norm and for $1\le p\le \infty$, the notation ${\parallel ·\parallel }_p$ stands for the Schatten p-norm defined as the p-norm of the vector of the singular values. Note that the Schatten ∞-norm, which is equal to the largest singular value coincides with the spectral norm $\parallel ·\parallel$; the Schatten 1-norm since it is the sum of the singular values is often referred to as the trace-norm; and the Schatten 2-norm coincides with the Frobenius norm. Schatten p-norms, as important special cases of unitarily invariant norms are treated in detail in a wonderful book by Bhatia [9].

Finally, the expression matrix-sequence refers to any sequence of the form ${\left\{A_n\right\}}_n$, where $A_n$ is a square matrix of size $d_n$ with $d_n$ strictly increasing so that $d_n\to \infty$ as $n\to \infty$. A r-block matrix-sequence, or simply a matrix-sequence if r can be deduced from context is a special ${\left\{A_n\right\}}_n$ in which the size of $A_n$ is $d_n=r{\phi }_n$, with $r\ge 1\in \mathbb{N}$ fixed and ${\phi }_n\in \mathbb{N}$ strictly increasing.

2.2. Multi-Index Notation

To effectively deal with multilevel structures it is necessary to use multi-indices, which are vectors of the form $\mathbf{i}=(i_1,\cdots ,i_d)\in {\mathbb{Z}}^d$. The related notation is listed below

• $\mathbf{0},\mathbf{1},\mathbf{2},\cdots$ are vectors of all zeroes, ones, twos, etc.
• $\mathbf{h}\le \mathbf{k}$ means that $h_r\le k_r$ for all $r=1,\cdots ,d$. In general, relations between multi-indices are evaluated componentwise.
• Operations between multi-indices, such as addition, subtraction, multiplication, and division, are also performed componentwise.
• The multi-index interval $[\mathbf{h},\mathbf{k}]$ is the set $\{\mathbf{j}\in {\mathbb{Z}}^d:\mathbf{h}\le \mathbf{j}\le \mathbf{k}\}$. We always assume that the elements in an interval $[\mathbf{h},\mathbf{k}]$ are ordered in the standard lexicographic manner

  $${[\cdots {\left[\left[(j_1,\cdots ,j_d)\right]j_d=h_d,\cdots ,k_d\right]}_{j_{d-1}=h_{d-1},\cdots ,k_{d-1}}\cdots ]}_{j_1=h_1,\cdots ,k_1}.$$
• $\mathbf{j}=\mathbf{h},\cdots ,\mathbf{k}$ means that $\mathbf{j}$ varies from $\mathbf{h}$ to $\mathbf{k}$, always following the lexicographic ordering.
• $\mathbf{m}\to \infty$ means that $min\left(\mathbf{m}\right)={min}_{j=1,\cdots ,d}{m}_j\to \infty$.
• The product of all the components of $\mathbf{m}$ is denoted as $\nu \left(\mathbf{m}\right):={\prod }_{j=1}^d{m}_j$.

A multilevel matrix-sequence is a matrix-sequence ${\left\{A_n\right\}}_n$ such that n varies in some infinite subset of $\mathbb{N}$, $\mathbf{n}=\mathbf{n}\left(n\right)$ is a multi-index in ${\mathbb{N}}^d$ depending on n, and $\mathbf{n}\to \infty$ when $n\to \infty$. This is typical of many approximations of differential operators in d dimensions.

2.3. Singular Value and Eigenvalue Distributions of a Matrix-Sequence

Let ${\mu }_k$ be the Lebesgue measure in ${\mathbb{R}}^k$. Throughout this work, all terminology from measure theory (such as “measurable set,” “measurable function,” “almost everywhere”, etc.) always refers to the Lebesgue measure. Let $C_c\left(\mathbb{R}\right)$ (resp., $C_c\left(\mathbb{C}\right)$) be the space of continuous complex-valued functions with bounded support defined on $\mathbb{R}$ (resp., $\mathbb{C}$). If $A\in {\mathbb{C}}^{n\times n}$, the singular values and eigenvalues of A are denoted by ${\sigma }_1\left(A\right),\cdots ,{\sigma }_n\left(A\right)$ and ${\lambda }_1\left(A\right),\cdots ,{\lambda }_n\left(A\right)$, respectively. The set of the eigenvalues (i.e., the spectrum) of A is denoted by $\Lambda \left(A\right)$.

Definition 1.

(Singular value and eigenvalue distribution of a matrix-sequence). Let ${\left\{A_n\right\}}_n$ be a matrix-sequence, with $A_n$ of size $d_n$, and let $\psi :D\subset {\mathbb{R}}^t\to {\mathbb{C}}^{r\times r}$ be a measurable function defined on a set D with $0<{\mu }_t\left(D\right)<\infty$.

• We say that ${\left\{A_n\right\}}_n$ has an (asymptotic) singular value distribution described by ψ, and we write ${\left\{A_n\right\}}_n{\sim }_{\sigma }\psi$, if

  $$\underset{n\to \infty }{lim}{\displaystyle \frac{1}{d_n}}\sum _{i=1}^{d_n}F\left({\sigma }_i\left(A_n\right)\right)={\displaystyle \frac{1}{{\mu }_t\left(D\right)}}{\int }_D{\displaystyle \frac{1}{r}}\sum _{i=1}^{r}F\left({\sigma }_i\left(\psi \left(\mathbf{x}\right)\right)\right)d\mathbf{x},\forall F\in C_c\left(\mathbb{R}\right).$$
• We say that ${\left\{A_n\right\}}_n$ has an (asymptotic) spectral (or eigenvalue) distribution described by ψ, and we write ${\left\{A_n\right\}}_n{\sim }_{\lambda }\psi$, if

  $$\underset{n\to \infty }{lim}{\displaystyle \frac{1}{d_n}}\sum _{i=1}^{d_n}F\left({\lambda }_i\left(A_n\right)\right)={\displaystyle \frac{1}{{\mu }_t\left(D\right)}}{\int }_D{\displaystyle \frac{1}{r}}\sum _{i=1}^{r}F\left({\lambda }_i\left(\psi \left(\mathbf{x}\right)\right)\right)d\mathbf{x},\forall F\in C_c\left(\mathbb{C}\right).$$
• If ψ describes both the singular value and eigenvalue distribution of ${\left\{A_n\right\}}_n$, we write ${\left\{A_n\right\}}_n{\sim }_{\sigma ,\lambda }\psi$.

In this case, the function ψ is referred to as the eigenvalue (or spectral) symbol of ${\left\{A_n\right\}}_n$.

The same definition when the considered matrix-sequence shows a multilevel structure. In that case n is replaced by $\mathbf{n}$ uniformly in $A_n$ and $d_n$.

The informal meaning behind the spectral distribution definition is as follows: if $\psi$ is continuous, then a suitable ordering of the eigenvalues ${\left\{{\lambda }_j\left(A_n\right)\right\}}_{j=1,\cdots ,d_n}$, assigned in correspondence with an equispaced grid on D, reconstructs approximately the r surfaces $\mathbf{x}\mapsto {\lambda }_i\left(\psi \left(\mathbf{x}\right)\right)$, $i=1,\cdots ,r$. For example, in the simplest case where $t=1$ and $D=[a,b]$, $d_n=nr$, the eigenvalues of $A_n$ are approximately equal — up to a few potential outliers — to ${\lambda }_i\left(\psi \left(x_j\right)\right)$, where

$$x_j=a+j{\displaystyle \frac{(b-a)}{n}},j=1,\cdots ,n,i=1,\cdots ,r.$$

If $t=2$ and $D=[a_1,b_1]\times [a_2,b_2]$, $d_n=n^{2}r$, the eigenvalues of $A_n$ are approximately equal — again up to a few potential outliers — to ${\lambda }_i\left(\psi (x_{j_1},y_{j_2})\right)$, where

$$x_{j_1}=a_1+j_1{\displaystyle \frac{(b_1-a_1)}{n}},y_{j_2}=a_2+j_2{\displaystyle \frac{(b_2-a_2)}{n}},j_1,j_2=1,\cdots ,n,i=1,\cdots ,r.$$

If the considered structure is two-level then the subscript is $\mathbf{n}=(n_1,n_2)$ and $d_n=n_1{n}_{2}r$.

Furthermore, for $t\ge 3$, a similar reasoning applies.

Finally we report an observation which is useful in the following derivations.

Remark 1.

The relation ${\left\{A_n\right\}}_n{\sim }_{\lambda }f$ and $\Lambda \left(A_n\right)\subseteq S$ for all n imply that the range of f is a subset of the closure $\overline{S}$ of S. In particular ${\left\{A_n\right\}}_n{\sim }_{\lambda }f$ and $A_n$ positive definite for all n imply that f is nonnegative definite almost everywhere, simply nonnegative almost everywhere if $r=1$. The same applies when a multilevel matrix-sequence ${\left\{A_{\mathbf{n}}\right\}}_{\mathbf{n}}$ is considered and similar statements hold when singular values are taken into account.

2.4. Approximating Classes of Sequences

In this subsection, we present the notion of the approximating class of sequences and a related key result.

Definition 2 

(Approximating class of sequences). Let ${\left\{A_n\right\}}_n$ be a matrix-sequence and let ${\left\{{\left\{B_{n,j}\right\}}_n\right\}}_j$ be a class of matrix-sequences, with $A_n$ and $B_{n,j}$ of size $d_n$. We say that ${\left\{{\left\{B_{n,j}\right\}}_n\right\}}_j$ is an approximating class of sequences (a.c.s.) for ${\left\{A_n\right\}}_n$ if the following condition is met: for every j, there exists $n_j$ such that, for every $n\ge n_j$,

$$A_n=B_{n,j}+R_{n,j}+N_{n,j},$$

$$\mathrm{rank}\left(R_{n,j}\right)\le c\left(j\right)d_n\mathrm{and}\parallel N_{n,j}\parallel \le \omega \left(j\right),$$

where $n_j$, $c\left(j\right)$, and $\omega \left(j\right)$ depend only on j, and

$$\underset{j\to \infty }{lim}c\left(j\right)=\underset{j\to \infty }{lim}\omega \left(j\right)=0.$$

${\left\{{\left\{B_{n,j}\right\}}_n\right\}}_j\stackrel{\mathrm{a}.\mathrm{c}.\mathrm{s}.\mathrm{wrt}j}{\to }{\left\{A_n\right\}}_n$ denotes that ${\left\{{\left\{B_{n,j}\right\}}_n\right\}}_j$ is an a.c.s. for ${\left\{A_n\right\}}_n$.

The following theorem represents the expression of a related convergence theory and it is a powerful tool used, for example, in the construction of the GLT *-Algebra

Theorem 1. 

Let ${\left\{A_n\right\}}_n$, ${\left\{B_{n,j}\right\}}_n$, with $j,n\in \mathbb{N}$, be matrix-sequences and let $\psi ,{\psi }_j:D\subset {\mathbb{R}}^d\to {\mathbb{C}}^{r\times r}$ be measurable functions defined on a set D with positive and finite Lebesgue measure. Suppose that:

1. 

${\left\{B_{n,j}\right\}}_n{\sim }_{\sigma }{\psi }_j$ for every j;

2. 

${\left\{{\left\{B_{n,j}\right\}}_n\right\}}_j\stackrel{\mathrm{a}.\mathrm{c}.\mathrm{s}.\mathrm{wrt}j}{\to }{\left\{A_n\right\}}_n$;

3. 

${\psi }_j\to \psi$ in measure.

Then

$${\left\{A_n\right\}}_n{\sim }_{\sigma }\psi .$$

Moreover, if all the involved matrices are Hermitian, the first assumption is replaced by ${\left\{B_{n,j}\right\}}_n{\sim }_{\lambda }{\psi }_j$ for every j, and the other two are left unchanged, then ${\left\{A_n\right\}}_n{\sim }_{\lambda }\psi .$

We end this section by observing that the same definition can be given and corresponding results (with obvious changes) hold, when the involved matrix-sequences show a multilevel structure. In that case n is replaced by $\mathbf{n}$ uniformly in $A_n$, $B_{n,j}$, $d_n$.

2.5. Matrix-Sequences with Explicit or Hidden (Asymptotic) Structure

In this subsection, we introduce the three types of matrix structures that constitute the basic building blocks of the GLT *-algebras. To be more specific, for any $d,r$ positive integers we consider the set of d-level r-block GLT matrix-sequences. For any $d,r$ positive integers, the considered set forms a *-algebra of matrix-sequences, which is maximal and is isometrically equivalent to the maximal *-algebra of $2d$-variate $r\times r$ matrix-valued measurable functions (with respect to the Lebesgue measure) defined canonically over ${[0,1]}^d\times {[-\pi ,\pi ]}^d$; see [4,7,8,19,20,21] and references therein.

The reduced version is essential when dealing with approximations of integro-differential operators (also in fractional versions) defined over general (non-Cartesian) domains. The idea was presented in [35,36] and it was exhaustively developed in [5], where the GLT symbols are again measurable functions defined over $\Omega \times {[-\pi ,\pi ]}^d$, with $\Omega$ Peano-Jordan measurable and contained in ${[0,1]}^d$. Als the reduced versions form maximal *-algebras isometrically equivalent to the corresponding maximal *-algebras of measurable functions. The considered GLT *-algebras represent rich examples of hidden (asymptotic) structures. Their building blocks are formed by two classes of explicit algebraic structures, d-level r-block Toeplitz and sampling diagonal matrix-sequences (see Section 2.7 and Section 2.8), plus the asymptotic structures given by the zero-distributed matrix-sequences; see Section 2.6. It is worth noticing that the latter class plays the role of compact operators with respect to bounded linear operators and in fact they form a two-sided ideal of matrix-sequences with respect to any of the GLT *-algebras.

2.6. Zero-Distributed Sequences

Zero-distributed sequences are defined as matrix-sequences ${\left\{A_n\right\}}_n$ such that ${\left\{A_n\right\}}_n{\sim }_{\sigma }0$. Note that, for any $r\ge 1$, ${\left\{A_n\right\}}_n{\sim }_{\sigma }0$ is equivalent to ${\left\{A_n\right\}}_n{\sim }_{\sigma }{O}_r$, where $O_r$ is the $r\times r$ zero matrix. The following theorem (see [34] and [20]), provides a useful characterization for detecting this type of sequence. We use the natural convention $1/\infty =0$.

Theorem 2. 

Let ${\left\{A_n\right\}}_n$ be a matrix-sequence, with $A_n$ of size $d_n$. Then

• ${\left\{A_n\right\}}_n{\sim }_{\sigma }0$ if and only if $A_n=R_n+N_n$ with $\mathrm{rank}\left(R_n\right)/d_n\to 0$ and $\parallel N_n\parallel \to 0$ as $n\to \infty$;
• ${\left\{A_n\right\}}_n{\sim }_{\sigma }0$ if there exists $p\in [1,\infty ]$ such that $\parallel A_n{\parallel }_p/d_n^{1/p}\to 0$ as $n\to \infty$.

As in Section 2.4, the same definition can be given and corresponding result (with obvious changes) holds, when the involved matrix-sequences show a multilevel structure. In that case n is replaced by $\mathbf{n}$ uniformly in $A_n,N_n,R_n,d_n$.

2.7. Multilevel Block Toeplitz Matrices

Given $\mathbf{n}\in {\mathbb{N}}^d$, a matrix of the form

$${\left[A_{\mathbf{i}-\mathbf{j}}\right]}_{\mathbf{i},\mathbf{j}=1}^{\mathbf{n}}\in {\mathbb{C}}^{\nu \left(\mathbf{n}\right)r\times \nu \left(\mathbf{n}\right)r},$$

with blocks $A_{\mathbf{k}}\in {\mathbb{C}}^{r\times r}$, $\mathbf{k}\in [-(\mathbf{n}-\mathbf{1}),\cdots ,\mathbf{n}-\mathbf{1}]$, is called multilevel block Toeplitz, or more precisely, d-level r-block Toeplitz matrix.

Given a matrix-valued function $f:{[-\pi ,\pi ]}^d\to {\mathbb{C}}^{r\times r}$ belonging to $L^1\left({[-\pi ,\pi ]}^d\right)$, the $\mathbf{n}$-th Toeplitz matrix associated with f is defined as

$$T_{\mathbf{n}}\left(f\right):={\left[{\widehat{f}}_{\mathbf{i}-\mathbf{j}}\right]}_{\mathbf{i},\mathbf{j}=\mathbf{1}}^{\mathbf{n}}\in {\mathbb{C}}^{\nu \left(\mathbf{n}\right)r\times \nu \left(\mathbf{n}\right)r},$$

where

$${\widehat{f}}_{\mathbf{k}}={\displaystyle \frac{1}{{\left(2\pi \right)}^d}}{\int }_{{[-\pi ,\pi ]}^d}f\left(\mathbf{\theta }\right)e^{-i(\mathbf{k},\mathbf{\theta })}d\mathbf{\theta }\in {\mathbb{C}}^{r\times r},\mathbf{k}\in {\mathbb{Z}}^d,$$

are the Fourier coefficients of f, in which i denotes the imaginary unit, the integrals are computed componentwise, and $\left(\mathbf{k},\mathbf{\theta }\right)=k_1{\theta }_1+\cdots +k_d{\theta }_d$. Equivalently, $T_{\mathbf{n}}\left(f\right)$ can be expressed as

$$T_{\mathbf{n}}\left(f\right)=\sum _{|j_1|<n_1}\cdots \sum _{|j_d|<n_d}{J}_{n_1}^{\left(j_1\right)}\otimes \cdots \otimes J_{n_d}^{\left(j_d\right)}\otimes {\widehat{f}}_{(j_1,\cdots ,j_d)},$$

where ⊗ denotes the Kronecker tensor product between matrices and $J_m^{\left(l\right)}$ is the matrix of order m whose $(i,j)$ entry equals 1 if $i-j=l$ and zero otherwise.

${\left\{T_{\mathbf{n}}\left(f\right)\right\}}_{\mathbf{n}\in {\mathbb{N}}^d}$ is the family of (multilevel block) Toeplitz matrices associated with f, which is called the generating function.

2.8. Block Diagonal Sampling Matrices

Given $d\ge 1$, $\mathbf{n}\in {\mathbb{N}}^d$, and a function $a:{[0,1]}^d\to {\mathbb{C}}^{r\times r}$, we define the multilevel block diagonal sampling matrix $D_{\mathbf{n}}\left(a\right)$ as the block diagonal matrix

$$D_{\mathbf{n}}\left(a\right)={\mathrm{diag}}_{\mathbf{i}=\mathbf{1},\cdots ,\mathbf{n}}a\left({\displaystyle \frac{\mathbf{i}}{\mathbf{n}}}\right)\in {\mathbb{C}}^{\nu \left(\mathbf{n}\right)r\times \nu \left(\mathbf{n}\right)r}.$$

2.9. The *-Algebra of d-Level r-Block GLT Matrix-Sequences

Let $r\ge 1$ be a fixed integer. A multilevel r-block GLT sequence, or simply a GLT sequence if we do not need to specify r, is a special multilevel r-block matrix-sequence equipped with a measurable function $\kappa :{[0,1]}^d\times {[-\pi ,\pi ]}^d\to {\mathbb{C}}^{r\times r}$, $d\ge 1$, called symbol. The symbol is essentially unique, in the sense that if $\kappa ,\varsigma$ are two symbols of the same GLT sequence, then $\kappa =\varsigma$ almost everywhere. We write ${\left\{A_n\right\}}_n{\sim }_{\mathrm{GLT}}\kappa$ to denote that ${\left\{A_n\right\}}_n$ is a GLT sequence with symbol $\kappa$.

It can be proven that the set of multilevel block GLT sequences is the *-algebra generated by the three classes of sequences defined in Section 2.6, Section 2.7, Section 2.8: zero-distributed, multilevel block Toeplitz, and block diagonal sampling matrix-sequences. The GLT class satisfies several algebraic and topological properties that are treated in detail in [7,8,20,21]. Here, we focus on the main operative properties, listed below, that represent a complete characterization of GLT sequences, equivalent to the full constructive definition.

GLT Axioms

• GLT 1. If ${\left\{A_{\mathbf{n}}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}\kappa$ then ${\left\{A_{\mathbf{n}}\right\}}_{\mathbf{n}}{\sim }_{\sigma }\kappa$ in the sense of Definition 1, with $t=2d$ and $D={[0,1]}^d\times {[-\pi ,\pi ]}^d$. Moreover, if each $A_{\mathbf{n}}$ is Hermitian, then ${\left\{A_{\mathbf{n}}\right\}}_{\mathbf{n}}{\sim }_{\lambda }\kappa$, again in the sense of Definition 1 with $t=2d$.
• GLT 2. We have

  -

  ${\left\{T_{\mathbf{n}}\left(f\right)\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}\kappa (\mathbf{x},\mathbf{\theta })=f(\mathbf{\theta })$ if $f:{[-\pi ,\pi ]}^d\to {\mathbb{C}}^{r\times r}$ is in $L^1\left({[-\pi ,\pi ]}^d\right)$;

  -

  ${\left\{D_{\mathbf{n}}\left(a\right)\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}\kappa (\mathbf{x},\mathbf{\theta })=a\left(\mathbf{x}\right)$ if $a:{[0,1]}^d\to {\mathbb{C}}^{r\times r}$ is Riemann-integrable;

  -

  ${\left\{Z_{\mathbf{n}}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}\kappa (\mathbf{x},\mathbf{\theta })=O_r$ if and only if ${\left\{Z_{\mathbf{n}}\right\}}_{\mathbf{n}}{\sim }_{\sigma }0$.

• GLT 3. If ${\left\{A_{\mathbf{n}}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}\kappa$ and ${\left\{B_{\mathbf{n}}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}\varsigma$, then:

  -

  ${\left\{A_{\mathbf{n}}^{*}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}{\kappa }^{*}$;

  -

  ${\{\alpha A_{\mathbf{n}}+\beta B_{\mathbf{n}}\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}\alpha \kappa +\beta \varsigma$ for all $\alpha ,\beta \in \mathbb{C}$;

  -

  ${\left\{A_{\mathbf{n}}{B}_{\mathbf{n}}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}\kappa \varsigma$;

  -

  ${\left\{A_{\mathbf{n}}^{†}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}{\kappa }^{-1}$, provided that $\kappa$ is invertible almost everywhere.

• GLT 4.${\left\{A_{\mathbf{n}}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}\kappa$ if and only if there exist ${\left\{B_{\mathbf{n},j}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}{\kappa }_j$ such that ${\left\{{\left\{B_{\mathbf{n},j}\right\}}_{\mathbf{n}}\right\}}_j\stackrel{\mathrm{a}.\mathrm{c}.\mathrm{s}.\mathrm{wrt}j}{\to }{\left\{A_{\mathbf{n}}\right\}}_{\mathbf{n}}$ and ${\kappa }_j\to \kappa$ in measure.
• GLT 5. If ${\left\{A_{\mathbf{n}}\right\}}_n{\sim }_{\mathrm{GLT}}\kappa$ and $A_{\mathbf{n}}=X_{\mathbf{n}}+Y_{\mathbf{n}}$, where

  -

  every $X_{\mathbf{n}}$ is Hermitian,

  -

  $\left|\right|X_{\mathbf{n}}\left|\right|,\left|\right|Y_{\mathbf{n}}\left|\right|\le C$ for some constant C independent of n,

  -

  $N{\left(\mathbf{n}\right)}^{-1}{\parallel Y_{\mathbf{n}}\parallel }_1\to 0$,

  then ${\left\{A_{\mathbf{n}}\right\}}_n{\sim }_{\lambda }\kappa$.
• GLT 6. If ${\left\{A_{\mathbf{n}}\right\}}_n{\sim }_{\mathrm{GLT}}\kappa$ and each $A_{\mathbf{n}}$ is Hermitian, then ${\left\{f\left(A_{\mathbf{n}}\right)\right\}}_n{\sim }_{\mathrm{GLT}}f\left(\kappa \right)$ for every continuous function $f:\mathbb{C}\to \mathbb{C}$.

Note that, by GLT 1, it is always possible to obtain the singular value distribution from the GLT symbol, while the eigenvalue distribution can only be deduced either if the involved matrices are Hermitian or the related matrix-sequence is quasi-Hermitian in the sense of GLT 5.

3. Geometric Mean of GLT Matrix-Sequences

This section discusses the geometric mean of positive definite matrices, starting with the well-established case of two matrices and then considering multiple matrices. In particular we give distribution results in the case of GLT matrix-sequences.

3.1. Means of Two Matrices

The geometric mean of two positive numbers a and b is simply $\sqrt{ab}$, a fact well-known from basic arithmetic. However, extending this concept to HPD matrices introduces a number of challenges, as matrix multiplication is not commutative. The question of how to define a geometric mean for matrices in a way that preserves key properties such as congruence invariance, consistency with scalars, and symmetry was solved by ALM [1]. Their work presented an axiomatic approach to defining the geometric mean of two HPD matrices previously defined in Equation (1). ALM formalized the geometric mean for matrices by establishing a set of ten essential properties, known as the ALM axioms. These axioms ensure the geometric mean behaves appropriately in the matrix setting. Here are three key properties:

1. Permutation invariance: $G(A,B)=G(B,A)$ for all $A,B\in {\mathcal{P}}_n$.

2. Congruence invariance: $G(M^{*}AM,M^{*}BM)=M^{*}G(A,B)M$ for all $A,B\in {\mathcal{P}}_n$ and all invertible matrices $M\in {\mathbb{C}}^{n\times n}$.

3. Consistency with scalars: $G(A,B)={\left(AB\right)}^{1/2}$ for all commuting $A,B\in {\mathcal{P}}_n$ (note that $AB\in {\mathcal{P}}_n$ for all commuting $A,B\in {\mathcal{P}}_n$ because ${\left(AB\right)}^{*}=B^{*}{A}^{*}=BA=AB$ and $AB$ is similar to the HPD matrix $A^{1/2}B{A}^{1/2}=A^{-1/2}\left(AB\right)A^{1/2}$).

When considering sequences of matrices, particularly in the framework of GLT sequences, the geometric mean operation is well-preserved under the structure of GLT sequences. If $r=d=1$ and we consider two scalar unilevel GLT matrix-sequences, that is, ${\left\{A_n\right\}}_n{\sim }_{\mathrm{GLT}}\kappa$ and ${\left\{B_n\right\}}_n{\sim }_{\mathrm{GLT}}\xi$, where $A_n,B_n\in {\mathcal{P}}_n$ are HPD matrices for every n, the matrix-sequence of their geometric mean ${\left\{G(A_n,B_n)\right\}}_n$ also forms a scalar unilevel GLT matrix-sequence. The symbol of the resulting sequence is the geometric mean of the individual symbols $\kappa$ and $\xi$.

Theorem 3 

([20], Theorem 10.2). Let $r=d=1$. Suppose ${\left\{A_n\right\}}_n{\sim }_{\mathrm{GLT}}\kappa$ and ${\left\{B_n\right\}}_n{\sim }_{\mathrm{GLT}}\xi$, where $A_n,B_n\in {\mathcal{P}}_n$ for every n. Assume that at least one between κ and ξ is nonzero almost everywhere. Then

$${\left\{G(A_n,B_n)\right\}}_n{\sim }_{\mathrm{GLT}}{\left(\kappa \xi \right)}^{1/2}$$

(7)

and

$${\left\{G(A_n,B_n)\right\}}_n{\sim }_{\sigma ,\lambda }{\left(\kappa \xi \right)}^{1/2}.$$

(8)

The previous result is easily extended to the case of matrix-sequence resulting from the geometric mean of two block multilevel GLT matrix-sequences, thanks to powerful *-algebra structures of considered spaces described in Section 2.9. Indeed the following two generalizations of Theorem 3 hold.

Theorem 4. 

Let $r=1$ and $d\ge 1$. Suppose ${\left\{A_{\mathbf{n}}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}\kappa$ and ${\left\{B_{\mathbf{n}}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}\xi$, where $A_{\mathbf{n}},B_{\mathbf{n}}\in {\mathcal{P}}_{\nu \left(\mathbf{n}\right)}$ for every multi-index $\mathbf{n}$. Assume that at least one between κ and ξ is nonzero almost everywhere. Then

$${\left\{G(A_{\mathbf{n}},B_{\mathbf{n}})\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}{\left(\kappa \xi \right)}^{1/2}$$

(9)

and

$${\left\{G(A_{\mathbf{n}},B_{\mathbf{n}})\right\}}_{\mathbf{n}}{\sim }_{\sigma ,\lambda }{\left(\kappa \xi \right)}^{1/2}.$$

(10)

Proof. Since both $A_{\mathbf{n}}$ and $B_{\mathbf{n}}$ are positive definite for every multi-index $\mathbf{n}$, the matrix-sequence ${\left\{G(A_{\mathbf{n}},B_{\mathbf{n}})\right\}}_{\mathbf{n}}$ is well defined according to formula (1) since $A_{\mathbf{n}}^{-1/2}$ is well defined for every multi-index $\mathbf{n}$. According to the assumption, we start with the case where $\kappa$ is nonzero almost everywhere. Hence the matrix-sequence ${\left\{A_{\mathbf{n}}^{-1}\right\}}_{\mathbf{n}}$ is a GLT matrix-sequence with GLT symbol ${\kappa }^{-1}$ by Axiom GLT 3, part 3. Since the square root is continuous and well defined over positive definite matrices also the matrix-sequence ${\left\{A_{\mathbf{n}}^{-1/2}\right\}}_{\mathbf{n}}$ is a GLT matrix-sequence with GLT symbol ${\kappa }^{-1/2}$ by virtue of Axiom GLT 6.

Now using two times GLT 3, part 2, we infer that ${\left\{A_{\mathbf{n}}^{-1/2}{B_{\mathbf{n}}{A}_{\mathbf{n}}}^{-1/2}\right\}}_{\mathbf{n}}$ is a GLT matrix-sequence with GLT symbol ${\kappa }^{-1}\xi$, where $X_{\mathbf{n}}=A_{\mathbf{n}}^{-1/2}{B_{\mathbf{n}}{A}_{\mathbf{n}}}^{-1/2}$ is positive definite because of the Sylvester inertia law. Hence the square root of $X_{\mathbf{n}}$ is well defined and by exploiting again Axiom GLT 6 we deduce that ${\left\{{\left(A_{\mathbf{n}}^{-1/2}{B_{\mathbf{n}}{A}_{\mathbf{n}}}^{-1/2}\right)}^{1/2}\right\}}_{\mathbf{n}}$ is a GLT matrix-sequence with GLT symbol ${\left({\kappa }^{-1}\xi \right)}^{1/2}$. Finally, by exploiting Axiom GLT 3, part 3, we have ${\left\{A_{\mathbf{n}}^{1/2}\right\}}_{\mathbf{n}}$ is a GLT matrix-sequence with GLT symbol ${\kappa }^{1/2}$ and the application of GLT 3, part 2, two time leads to the desired conclusion

$${\left\{G(A_{\mathbf{n}},B_{\mathbf{n}})\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}{\left(\kappa \xi \right)}^{1/2},$$

where the latter and Axiom GLT 3 imply ${\left\{G(A_{\mathbf{n}},B_{\mathbf{n}})\right\}}_{\mathbf{n}}{\sim }_{\sigma ,\lambda }{\left(\kappa \xi \right)}^{1/2}$.

Finally the other case where $\xi$ is nonzero almost everywhere has the very same proof since $G(·,·)$ is invariant under permutations and hence

$$G(A,B)=A^{1/2}{\left(A^{-1/2}B{A}^{-1/2}\right)}^{1/2}{A}^{1/2}=G(B,A)=B^{1/2}{\left(B^{-1/2}A{B}^{-1/2}\right)}^{1/2}{B}^{1/2},$$

so that the the same steps can be repeated by exchanging $A_{\mathbf{n}}$ and $B_{\mathbf{n}}$. •

Theorem 5. 

Let $r,d\ge 1$. Suppose ${\left\{A_{\mathbf{n}}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}\kappa$ and ${\left\{B_{\mathbf{n}}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}\xi$, where $A_{\mathbf{n}},B_{\mathbf{n}}\in {\mathcal{P}}_{\nu \left(\mathbf{n}\right)}$ for every multi-index $\mathbf{n}$. Assume that at least one between the minimal eigenvalue of κ and the minimal eigenvalue of ξ is nonzero almost everywhere. Then

$${\left\{G(A_{\mathbf{n}},B_{\mathbf{n}})\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}G(\kappa ,\xi )$$

(11)

and

$${\left\{G(A_{\mathbf{n}},B_{\mathbf{n}})\right\}}_{\mathbf{n}}{\sim }_{\sigma ,\lambda }G(\kappa ,\xi ).$$

(12)

Furthermore $G(\kappa ,\xi )={\left(\kappa \xi \right)}^{1/2}$ whenever the GLT symbols κ and ξ commute.

Proof. The case of $r=1$ is already contained in Theorem 4, so we assume $r>1$ i.e. a true GLT block setting. The proof is in fact a repetition of that of the previous theorem with the only attention to GLT symbol part where the multiplication is noncommutative.

Since both $A_{\mathbf{n}}$ and $B_{\mathbf{n}}$ are positive definite for every multi-index $\mathbf{n}$, the matrix-sequence ${\left\{G(A_{\mathbf{n}},B_{\mathbf{n}})\right\}}_{\mathbf{n}}$ is well defined according to formula (1) since $A_{\mathbf{n}}^{-1/2}$ is well defined for every multi-index $\mathbf{n}$. According to the assumption, we start with the case where $\kappa$ is invertible almost everywhere so that ${\left\{A_{\mathbf{n}}^{-1}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}{\kappa }^{-1}$ by Axiom GLT 3, part 3, and ${\left\{A_{\mathbf{n}}^{-1/2}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}{\kappa }^{-1/2}$ thanks to Axiom GLT 6.

Now using two times GLT 3, part 2, we have ${\left\{A_{\mathbf{n}}^{-1/2}{B_{\mathbf{n}}{A}_{\mathbf{n}}}^{-1/2}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}{\kappa }^{-1/2}\xi {\kappa }^{-1/2}$, where

$$X_{\mathbf{n}}=A_{\mathbf{n}}^{-1/2}{B_{\mathbf{n}}{A}_{\mathbf{n}}}^{-1/2}$$

is positive definite because of the Sylvester inertia law. Hence the square root of $X_{\mathbf{n}}$ is well defined and by exploiting again Axiom GLT 6 we obtain ${\left\{{\left(A_{\mathbf{n}}^{-1/2}{B_{\mathbf{n}}{A}_{\mathbf{n}}}^{-1/2}\right)}^{1/2}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}{\left({\kappa }^{-1/2}\xi {\kappa }^{-1/2}\right)}^{1/2}$. Finally, by exploiting Axiom GLT 3, part 3, on the matrix-sequence ${\left\{A_{\mathbf{n}}^{1/2}\right\}}_{\mathbf{n}}$ and using GLT 3, part 2, two times, we conclude

$${\left\{G(A_{\mathbf{n}},B_{\mathbf{n}})\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}{\kappa }^{1/2}{\left({\kappa }^{-1/2}\xi {\kappa }^{-1/2}\right)}^{1/2}{\kappa }^{1/2},$$

where the symbol ${\kappa }^{1/2}{\left({\kappa }^{-1/2}\xi {\kappa }^{-1/2}\right)}^{1/2}{\kappa }^{1/2}$ is exactly $G(\kappa ,\xi )$. Furthermore, relation (11) and Axiom GLT 3 imply ${\left\{G(A_{\mathbf{n}},B_{\mathbf{n}})\right\}}_{\mathbf{n}}{\sim }_{\sigma ,\lambda }G(\kappa ,\xi )$ where $G(\kappa ,\xi )={\left(\kappa \xi \right)}^{1/2}$, whenever $\kappa$ and $\xi$ commute.

Finally the remaining case where $\xi$ is invertible almost everywhere has the very same proof since $G(·,·)$ is invariant under permutations and hence

$$G(A,B)=A^{1/2}{\left(A^{-1/2}B{A}^{-1/2}\right)}^{1/2}{A}^{1/2}=G(B,A)=B^{1/2}{\left(B^{-1/2}A{B}^{-1/2}\right)}^{1/2}{B}^{1/2},$$

so that the the same steps can be repeated by exchanging $A_{\mathbf{n}}$ and $B_{\mathbf{n}}$ and by exchanging $\kappa$ and $\xi$. •

3.2. Mean of More Than Two Matrices

In this section, we describe the iterative method used to compute the Karcher mean for more than two HPD matrices. The Karcher mean is an extension of the geometric mean to more than two matrices and can be computed using an iterative method based on the Richardson-like iteration. As detailed earlier in the introduction (see Equation 6), the iteration updates the approximation at each step based on the logarithmic correction term. The step-size parameter $\theta$ is dynamically computed during the iteration, ensuring that the process converges efficiently by accounting for the condition numbers of the matrices involved.

Specifically, $\theta$ is given by

$$\theta ={\displaystyle \frac{2}{\gamma +\beta }},$$

(13)

with $\beta$ and $\gamma$ computed as

$$[\beta ,\gamma ]=\left[\sum _{j=1}^k{\displaystyle \frac{log\left(c_j\right)}{c_j-1}},\sum _{j=1}^k{c}_j{\displaystyle \frac{log\left(c_j\right)}{c_j-1}}\right],$$

(14)

where $c_j={\displaystyle \frac{{\lambda }_{\max }\left(M_j\right)}{{\lambda }_{\min }\left(M_j\right)}},$ and $M_j=G^{1/2}{A}_j^{-1}{G}^{1/2},$G is the current approximation of the Karcher mean.

The Richardson-like iteration can be implemented in different equivalent forms

$$\begin{array}{cc}\hfill X& =X-\theta X^{1/2}\left(\sum _{i=1}^{k}log\left(X^{1/2}{A}_i^{-1}{X}^{1/2}\right)\right)X^{1/2},G=X,\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill X& =X^{1/2}\left(I-\theta \sum _{i=1}^{k}log\left(X^{1/2}{A}_i^{-1}{X}^{1/2}\right)\right)X^{1/2},G=X,\hfill \end{array}$$

(16)

$$\begin{array}{cc}\hfill Y& =Y-\theta Y^{1/2}\left(\sum _{i=1}^{k}log\left(Y^{1/2}{A}_i{Y}^{1/2}\right)\right)Y^{1/2},G=Y^{-1},\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill Y& =Y^{1/2}\left(I-\theta \sum _{i=1}^{k}log\left(Y^{1/2}{A}_i{Y}^{1/2}\right)\right)Y^{1/2},G=Y^{-1}.\hfill \end{array}$$

(18)

Among these equivalent formulations, the first one, Equation (15), is the most practical for implementation. It avoids matrix inversions, which can introduce numerical instabilities and increase computational complexity. While formulations Equation () and Equation () aim to reduce the number of matrix inversions, the final step requires inverting the result, which can lead to inaccuracies, especially for poorly conditioned matrices. Additionally, Equation () retains the simplicity of direct matrix operations without introducing unnecessary complications. A numerically more efficient approach uses Cholesky factorization to reduce the computational cost of forming matrix square roots at every step, enhancing efficiency as forming the Cholesky factor costs less than computing a full matrix square root [24].

Suppose $X_{\nu }=R_{\nu }^T{R}_{\nu }$ is the Cholesky decomposition of $X_{\nu }$, where $R_{\nu }$ is an upper triangular matrix. The iteration step can be rewritten as

$$X_{\nu +1}=X_{\nu }+\theta R_{\nu }^T\left(\sum _{i=1}^{k}log\left(R_{\nu }^{-T}{A}_i{R}_{\nu }^{-1}\right)\right)R_{\nu }.$$

(19)

In this formulation, the Cholesky factor $R_{\nu }$ is updated at each iteration. The condition number of the Cholesky factor $R_{\nu }$ in the spectral norm is the square root of the condition number of $X_{\nu }$, thus ensuring good numerical accuracy. For this heuristic to be effective, it is essential that $X_0$ provides a good approximation of G. Therefore, selecting $X_0$ as the Cheap mean is critical. An adaptive version of this iteration has been proposed and implemented in the Matrix Means Toolbox [11].

Of course the Richardson-like iteration is relevant for computing efficiently the Karcher mean: we also exploit its formal expression for theoretical purposes when dealing with sequences of matrices, particularly those involving GLT sequences. For the theory we come back at relation (15) and we consider the associated iteration

$$X_{\nu +1}=X_{\nu }-\theta X_{\nu }^{1/2}\left(\sum _{i=1}^{k}log\left(X_{\nu }^{1/2}{A}_i^{-1}{X}_{\nu }^{1/2}\right)\right)X_{\nu }^{1/2},$$

(20)

with $X_0$ given positive definite matrix. We know that $X_{\nu }$ converges to the geometric mean of $A_1,\dots ,A_k$ as $\nu$ tends to infinity for every fixed positive definite initial guess $X_0$.

Fix $r,d\ge 1$. Suppose now that the block multitivel sequence of matrices ${\left\{A_{\mathbf{n}}^{\left(i\right)}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}{\kappa }_i$ for $i=1,\cdots ,k$ are given, where $A_{\mathbf{n}}^{\left(1\right)},\cdots ,A_n^{\left(k\right)}$ are positive definite for every multi-index $\mathbf{n}$.

Due to the positive definiteness of the matrices ${\left\{A_{\mathbf{n}}^{\left(i\right)}\right\}}_{\mathbf{n}}$ and because of ${\left\{A_{\mathbf{n}}^{\left(i\right)}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}{\kappa }_i$, from Axiom GLT 1 it follows that each ${\kappa }_i$ is nonnegative definite almost everywhere (see Remark 1).

In this setting, it is conjectured that the sequence of Karcher means ${\left\{G(A_{\mathbf{n}}^{\left(1\right)},\cdots ,A_{\mathbf{n}}^{\left(k\right)})\right\}}_{\mathbf{n}}$ forms a new GLT matrix-sequence whose symbol is the geometric mean of the individual symbols ${\kappa }_1,\cdots ,{\kappa }_k$, specifically ${({\kappa }_1\cdots {\kappa }_k)}^{1/k}$ if all symbols commute and $G({\kappa }_1,\cdots ,{\kappa }_k)$ in the general case. In order to attack the problem the initial guess matrix-sequence ${\left\{X_{\mathbf{n},0}\right\}}_{\mathbf{n}}$ must be of GLT type with nonnegative definite GLT symbol. In this way thanks to (20) and using the GLT axioms in the way it is done in Theorem 5, we deduce easily that ${\left\{X_{\mathbf{n},\nu }\right\}}_{\mathbf{n}}$ is still a GLT matrix-sequence with symbol $g_{\nu }$ converging to $G({\kappa }_1,\cdots ,{\kappa }_k)$.

Finally Theorem 1 and Axiom GLT 4 could be applied if we prove that ${\left\{{\left\{X_{\mathbf{n},\nu }\right\}}_{\mathbf{n}}\right\}}_{\nu }$ if an a.c.s. for the limit sequence ${\left\{G(A_{\mathbf{n}}^{\left(1\right)},\cdots ,A_{\mathbf{n}}^{\left(k\right)})\right\}}_{\mathbf{n}}$. This could be proven using Schatten estimates like those in the second item of Theorem 2, but at the moment this is not easy because the known convergence proofs for the Karcher iterations are all based on pointwise convergence.

4. Numerical Experiments

In this section, we present and critically analyze several selected examples, by considering matrix-sequences of geometric means k HPD matrices for $k=2,3$. In particular our numerics show asymptotical spectral properties of the resulting matrix-sequences in accordance with the theoretical results (and conjectures) of the previous section. We introduce few examples in which the input matrix-sequences are either of Toeplitz type or are general r-block d-level GLT matrix-sequences, stemming from the approximations of differential operators via local methods such finite differences, finite elements, isogeometric analysis: in the first group we explore the geometric mean of two matrix-sequences and in the second group we consider the Karcher mean of three matrix-sequences, both taking into account one-dimensional (1D) and two-dimensional (2D) settings, i.e., $d=1,2$ and $r=1$; in the final group we deal with r-block GLT matrix-sequences with $r=2$. We anticipate the strong agreement of the numerical evidences with the theoretical results in Theorem 3, Theorem 4, Theorem 5, and with the conjecture regarding the Karcher mean, also when the matrix-sizes are quite moderate. The latter is a nontrivial numerical finding since all the theoretical results have an asymptotic spectral nature.

4.1. Example 1 (1D)

Let $A_n=T_n\left({(2-2cos\left(\theta \right))}^2\right)$ according to Section 2.7 and let $B_n$ be the finite difference discretization of the differential operator ${\left(\alpha \left(x\right)u^{\prime \prime }\right)}^{\prime \prime }$ on the interval $(0,1)$, with boundary conditions $u\left(0\right)=u\left(1\right)=u^{\prime }\left(0\right)=u^{\prime }\left(1\right)=0$, where $\alpha \left(x\right)$ is positive on $(0,1)$. For the fourth order boundary value problem

$${\left(\alpha \right(x\left)u^{\prime \prime }\right)}^{\prime \prime }=f\left(x\right),u\left(0\right)=u\left(1\right)=u^{\prime }\left(0\right)=u^{\prime }\left(1\right)=0,$$

we approximate the derivative $u^{\left(4\right)}\left(x\right)$ by using the second-order central FD scheme characterized by the stencil $(1,-4,6,-4,1)$. More specifically, for $\alpha ,u$ smooth enough, we have ${\left.{\left(\alpha \left(x\right)u^{\prime \prime }\left(x\right)\right)}^{\prime \prime }\right|}_{x=x_i}$ equal to

$${\displaystyle \frac{{\alpha }_{i-2}{u}_{i-2}-2({\alpha }_{i-1}+{\alpha }_i)u_{i-1}+({\alpha }_{i-1}+4{\alpha }_i+{\alpha }_{i+1})u_i-2({\alpha }_{i+1}+{\alpha }_i)u_{i+1}+{\alpha }_{i+1}{u}_{i+2}}{h^4}}+O\left(h^2\right),$$

for all $i=2,3,....,n+1$; here, $h={\displaystyle \frac{1}{n+3}}$ and $x_i=ih$ for $i=0,1,....,n+3$.

By taking into account the homogeneous boundary conditions and by neglecting the $O\left(h^2\right)$ approximation error, we approximate the nodal value $u\left(x_i\right)$ with the value of $u_i$ for $i=0,1,....,n+3$, where $u_0=u_1=u_{n+2}=u_{n+3}=0$ and $u={(u_2,....,u_{n+1})}^T$ is the solution of the linear system

$${\alpha }_{i-2}{u}_{i-2}-2({\alpha }_{i-1}+{\alpha }_i)u_{i-1}+({\alpha }_{i-1}+4{\alpha }_i+{\alpha }_{i+1})u_i-2({\alpha }_{i+1}+{\alpha }_i)u_{i+1}+{\alpha }_{i+1}{u}_{i+2}=h^{4}f\left(x\right),$$

for all $i=2,3,....,n+1$.

The structure of the resulting matrix $B_n=B_n\left(\alpha \right)$ is as reported below

(21)

Looking at $B_n$ in (21), we observe that it can be written as

$$B_n=B_n\left(\alpha \right)=D_n^{+}{K}_n^{+}+D_n{K}_n+D_n^{-}{K}_n^{-},$$

with

$$K_n^{+}=\left(\begin{array}{ccccccc}1& -2& 1& & & & \\ & 1& -2& 1& & & \\ & & 1& -2& 1& & \\ & & & 1& -2& \ddots & \\ & & & & 1& \ddots & 1\\ & & & & & \ddots & -2\\ & & & & & & 1\end{array}\right)=T_n(1-2e^{-i\theta }+e^{-2i\theta }),$$

$$K_n=\left(\begin{array}{ccccccc}4& -2& & & & & \\ -2& 4& -2& & & & \\ & -2& 4& -2& & & \\ & & -2& 4& -2& & \\ & & & \ddots & \ddots & \ddots & \\ & & & & -2& 4& -2\\ & & & & & -2& 4\end{array}\right)=T_n(4-2e^{-i\theta }-2e^{i\theta }),$$

$$K_n^{-}=\left(\begin{array}{ccccccc}1& & & & & & \\ -2& 1& & & & & \\ 1& -2& 1& & & & \\ & 1& -2& 1& & & \\ & & \ddots & \ddots & \ddots & & \\ & & & 1& -2& 1& \\ & & & & 1& -2& 1\end{array}\right)=T_n(1-2e^{i\theta }+e^{2i\theta }).$$

$$\begin{array}{ccc}\hfill D_n^{+}& =& \underset{i=1,2....,n}{\mathrm{diag}}{\alpha }_{i+2}=\underset{i=1,2....,n}{\mathrm{diag}}\alpha \left(x_{i+2}\right),\hfill \\ \hfill D_n& =& \underset{i=1,2....,n}{\mathrm{diag}}{\alpha }_{i+1}=\underset{i=1,2....,n}{\mathrm{diag}}\alpha \left(x_{i+1}\right),\hfill \\ \hfill D_n^{-}& =& \underset{i=1,2....,n}{\mathrm{diag}}{\alpha }_i=\underset{i=1,2....,n}{\mathrm{diag}}\alpha \left(x_i\right).\hfill \end{array}$$

It is easy to check that grids ${\mathcal{G}}_n^{+}={\left\{x_{i+2}\right\}}_{i=1,....,n},{\mathcal{G}}_n={\left\{x_{i+1}\right\}}_{i=1,....,n}$ and ${\mathcal{G}}_n^{-}={\left\{x_i\right\}}_{i=1,....,n}$ are asymptotically uniform in $[0,1]$, according to the notion given in [6]. In fact, for ${\mathcal{G}}_n^{+}={\left\{x_{i+2}\right\}}_{i=1,....,n}$ we have

$$\underset{i=1,2,\cdots ,n}{\max }\left|x_{i,n}-{\displaystyle \frac{i}{n}}\right|=\underset{i=1,2,\cdots ,n}{\max }\left|x_{i+2}-{\displaystyle \frac{i}{n}}\right|=\underset{i=1,2,\cdots ,n}{\max }\left|(i+2)h-{\displaystyle \frac{i}{n}}\right|.$$

$$.=\underset{i=1,2,\cdots ,n}{\max }\left|{\displaystyle \frac{i+2}{n+3}}-{\displaystyle \frac{i}{n}}\right|=\underset{i=1,2,\cdots ,n}{\max }\left|{\displaystyle \frac{ni+2n-i(n+3)}{n(n+3)}}\right|=\underset{i=1,2,\cdots ,n}{\max }\left|{\displaystyle \frac{-i}{n(n+3)}}\right|\le {\displaystyle \frac{n}{n(n+3)}}\to 0.$$

For ${\mathcal{G}}_n={\left\{x_{i+1}\right\}}_{i=1,....,n}$ we have

$$\underset{i=1,2,\cdots ,n}{\max }\left|x_{i,n}-{\displaystyle \frac{i}{n}}\right|=\underset{i=1,2,\cdots ,n}{\max }\left|x_{i+1}-{\displaystyle \frac{i}{n}}\right|=\underset{i=1,2,\cdots ,n}{\max }\left|(i+1)h-{\displaystyle \frac{i}{n}}\right|.$$

$$.=\underset{i=1,2,\cdots ,n}{\max }\left|{\displaystyle \frac{i+1}{n+3}}-{\displaystyle \frac{i}{n}}\right|=\underset{i=1,2,\cdots ,n}{\max }\left|{\displaystyle \frac{ni+n-i(n+3)}{n(n+3)}}\right|=\underset{i=1,2,\cdots ,n}{\max }\left|{\displaystyle \frac{-2i}{n(n+3)}}\right|\le {\displaystyle \frac{2n}{n(n+3)}}\to 0.$$

For ${\mathcal{G}}_n^{-}={\left\{x_i\right\}}_{i=1,....,n}$ we have

$$\underset{i=1,2,\cdots ,n}{\max }\left|x_{i,n}-{\displaystyle \frac{i}{n}}\right|=\underset{i=1,2,\cdots ,n}{\max }\left|x_i-{\displaystyle \frac{i}{n}}\right|=\underset{i=1,2,\cdots ,n}{\max }\left|ih-{\displaystyle \frac{i}{n}}\right|.$$

$$.=\underset{i=1,2,\cdots ,n}{\max }\left|{\displaystyle \frac{i}{n+3}}-{\displaystyle \frac{i}{n}}\right|=\underset{i=1,2,\cdots ,n}{\max }\left|{\displaystyle \frac{ni-i(n+3)}{n(n+3)}}\right|=\underset{i=1,2,\cdots ,n}{\max }\left|{\displaystyle \frac{-3i}{n(n+3)}}\right|\le {\displaystyle \frac{3n}{n(n+3)}}\to 0.$$

Hence, by $\mathbf{GLT}\mathbf{2}$, part 2, and by [6], we deduce that ${\left\{D_n^{+}\right\}}_n{\sim }_{\mathrm{GLT}}\alpha \left(x\right)$, ${\left\{D_n\right\}}_n{\sim }_{\mathrm{GLT}}\alpha \left(x\right)$ and ${\left\{D_n^{-}\right\}}_n{\sim }_{\mathrm{GLT}}\alpha \left(x\right)$. In conclusion, by invoking $\mathbf{GLT}\mathbf{2}-\mathbf{GLT}\mathbf{3}$, we infer

$${\left\{B_n\right\}}_n{\sim }_{\mathrm{GLT}}\alpha \left(x\right)(1-2e^{-i\theta }+e^{-2i\theta })+\alpha \left(x\right)(4-2e^{-i\theta }-2e^{i\theta })+\alpha \left(x\right)(1-2e^{i\theta }+e^{2i\theta })=\alpha \left(x\right){(2-2cos\left(\theta \right))}^2.$$

Eigenvalue Distribution

We begin by numerically verifying the eigenvalue distribution of the matrix-sequence ${\left\{G(A_n,B_n)\right\}}_n$ with respect to its GLT symbol ${\left(\kappa \xi \right)}^{1/2}$ according to Theorem 8, with $A_n=T_n\left({(2-2cos\left(\theta \right))}^2\right)$, $B_n=B_n\left(\alpha \right)$ as in (21) with $\alpha \left(x\right)=x$, $\kappa (x,\theta )={(2-2cos\left(\theta \right))}^2$, ${\xi (x,\theta ))=\alpha \left(x\right)(2-2cos\left(\theta \right))}^2$. In Figure 1, we compare the eigenvalues of the geometric mean with a uniform sampling of the symbol. It is evident that, as n increases, the symbol provides a better and better approximation of the eigenvalues. Similar results for the two-dimensional case are shown in Example Section 4.2, Figure 2.

Figure 1. Comparison between the symbol ${\left(\kappa \xi \right)}^{{\displaystyle \frac{1}{2}}}$ and $\mathrm{eig}\left(G(A_n,B_n)\right)$.

Figure 1. Comparison between the symbol ${\left(\kappa \xi \right)}^{{\displaystyle \frac{1}{2}}}$ and $\mathrm{eig}\left(G(A_n,B_n)\right)$.

4.2. Example 2 (2D)

Let $A_{\mathbf{n}}=T_{\mathbf{n}}\left(F({\theta }_1,{\theta }_2)\right)$ according to the notation in Section 2.7 with $d=2$ and $B_{\mathbf{n}}$ be the finite difference discretization of the differential operator

$${\displaystyle \frac{{\partial }^2}{\partial x^2}}\left(a(x,y){\displaystyle \frac{{\partial }^2}{\partial x^2}}\right)+{\displaystyle \frac{{\partial }^2}{\partial y^2}}\left(b(x,y){\displaystyle \frac{{\partial }^2}{\partial y^2}}\right),$$

on the open domain $\Omega ={(0,1)}^2$, with $a(x,y),b(x,y)$ nonnegative on the closure of $\Omega$, with homogeneous Dirichlet boundary conditions on $u(x,y)$, and zero normal derivatives at $\partial \Omega$.

Regarding $A_{\mathbf{n}}$, the generating function $F({\theta }_1,{\theta }_2)$ is given by

$$F({\theta }_1,{\theta }_2)=f\left({\theta }_1\right)+f\left({\theta }_2\right),$$

with $f\left(\theta \right)={(2-2cos\left(\theta \right))}^2$.

Thus, we have

$$A_{\mathbf{n}}=T_{\mathbf{n}}(f\left({\theta }_1\right)+f\left({\theta }_2\right))=T_{\mathbf{n}}(8+4{cos}^2\left({\theta }_1\right)+4{cos}^2\left({\theta }_2\right)-8cos\left({\theta }_1\right)-8cos\left({\theta }_1\right)).$$

As in the one-dimensional setting in Example Section 4.1, we apply the second-order central finite difference scheme separately to the x- and y-directions, in perfect analogy with the 1D case, and take $a(x,y)=\alpha \left(x\right)$, $b(x,y)=\alpha \left(y\right)$. Choosing $\alpha \left(t\right)=t$ the related problem is semielliptic, but it is of separable nature. This separable nature is reflected algebraically in a tensor decomposition of the whole approximation and in fact we find that resulting global matrix is given by the Kronecker sum

$$B_{\mathbf{n}}=B_{(n_1,n_2)}=B_{n_1}\left(\alpha \right)\otimes I_{n_2}+I_{n_1}\otimes B_{n_2}\left(\alpha \right),$$

where $I_n$ denotes the identity matrix of size n and $B_n\left(\alpha \right)$ is exactly the structure displayed in (21). By exploiting the GLT analysis in the one-dimensional case, the symbol for the two-dimensional matrix-sequence can be derived similarly. More precisely, we have

$${\left\{B_{\mathbf{n}}\right\}}_{\mathbf{n}}{\sim }_{\mathrm{GLT}}\alpha \left(x\right){(2-2cos\left({\theta }_1\right))}^2+\alpha \left(y\right){(2-2cos\left({\theta }_2\right))}^2.$$

As shown in Figure 2, the agreement is remarkable even for the largest eigenvalues for which the absolute discrepancy is higher.

Figure 2. Comparison between the symbol ${\left(\kappa \xi \right)}^{{\displaystyle \frac{1}{2}}}$ and $\mathrm{eig}(G(A_{\mathbf{n}},B_{\mathbf{n}})$.

Figure 2. Comparison between the symbol ${\left(\kappa \xi \right)}^{{\displaystyle \frac{1}{2}}}$ and $\mathrm{eig}(G(A_{\mathbf{n}},B_{\mathbf{n}})$.

4.3. Example 3 (1D)

Let $A_n^{\left(1\right)}=T_n(3+2cos\left(\theta \right))$, $A_n^{\left(2\right)}=D_n\left(a\right)$, and $A_n^{\left(3\right)}=D_n\left(a\right)T_n(4-2cos\left(\theta \right))D_n\left(a\right)$, where $T_n(·)$ is the Toeplitz operator for $d=1$ as in Section 2.7 and $D_n\left(a\right)$ is the diagonal matrix generated by the continuous function $a\left(x\right)=x^2$, according to the notation in Section 2.8.

Eigenvalue distribution

First, we aim to verify the eigenvalue distribution of more than two matrices using the Karcher mean. In Figure 3, we compare the eigenvalues of the Karcher mean with a uniform sampling of the resulting limit GLT symbol. It can be observed that the symbol provides a reliable approximation of the eigenvalues, and in fact, as n increases, the spectral distribution holds asymptotically. Similar results for the two-dimensional case are shown in Example Section 4.4, Figure 4. Both types of result corroborate the conjecture on the GLT nature of the limit matrix-sequence of the Karcher means, as discussed at the end of Section 3.2.

Figure 3. Comparison between the symbol ${\left({\kappa }_1{\kappa }_2{\kappa }_3\right)}^{{\displaystyle \frac{1}{3}}}$ and $\mathrm{eig}\left(G(A_n^{\left(1\right)},A_n^{\left(2\right)},A_n^{\left(3\right)})\right)$.

Figure 3. Comparison between the symbol ${\left({\kappa }_1{\kappa }_2{\kappa }_3\right)}^{{\displaystyle \frac{1}{3}}}$ and $\mathrm{eig}\left(G(A_n^{\left(1\right)},A_n^{\left(2\right)},A_n^{\left(3\right)})\right)$.

Figure 4. Comparison between the symbol ${\left({\kappa }_1{\kappa }_2{\kappa }_3\right)}^{{\displaystyle \frac{1}{3}}}$ and $\mathrm{eig}\left(G(A_{\mathbf{n}}^{\left(1\right)},A_{\mathbf{n}}^{\left(2\right)},A_{\mathbf{n}}^{\left(3\right)})\right)$.

Figure 4. Comparison between the symbol ${\left({\kappa }_1{\kappa }_2{\kappa }_3\right)}^{{\displaystyle \frac{1}{3}}}$ and $\mathrm{eig}\left(G(A_{\mathbf{n}}^{\left(1\right)},A_{\mathbf{n}}^{\left(2\right)},A_{\mathbf{n}}^{\left(3\right)})\right)$.

4.4. Example 4 (2D)

Consider the matrix $A_{\mathbf{n}}^{\left(1\right)}=T_{\mathbf{n}}\left(F({\theta }_1,{\theta }_2)\right)$, where $T_{\mathbf{n}}(·)$ is the Toeplitz operator as in Section 2.7 with $d=2$, and the function $F({\theta }_1,{\theta }_2)$ is defined as $F({\theta }_1,{\theta }_2)=f\left({\theta }_1\right)+f\left({\theta }_2\right)$, with $f\left(\theta \right)=3+2cos\left(\theta \right)$, resulting in $A_{\mathbf{n}}^{\left(1\right)}=T_{\mathbf{n}}(6+2cos\left({\theta }_1\right)+2cos\left({\theta }_2\right))$.

Also, consider the diagonal matrix $A_{\mathbf{n}}^{\left(2\right)}=D_{\mathbf{n}}\left(a\right)$, where, according to the notation in Section 2.8, $D_{\mathbf{n}}\left(a\right)$ is the diagonal sampling matrix generated by a continuous function $a(x,y)=x^2+y^2$, which is positive on the domain ${(0,1)}^2$.

We take the matrix $A_{\mathbf{n}}^{\left(3\right)}=D_{\mathbf{n}}\left(b\right)T_{\mathbf{n}}\left(G({\theta }_1,{\theta }_2)\right)D_{\mathbf{n}}\left(b\right)$, where $b(x,y)=1/x+1/y$ is positive and unbounded on the domain ${(0,1)}^2$, and where the generating function $G({\theta }_1,{\theta }_2)=f\left({\theta }_1\right)+f\left({\theta }_2\right)$, with $f\left(\theta \right)=4-2cos\left(\theta \right)$ implies that $A_{\mathbf{n}}^{\left(3\right)}$ is $T_{\mathbf{n}}(8-2cos\left({\theta }_1\right)-2cos\left({\theta }_2\right))$.

Also in the current example the agreement between the limit GLT symbol and the displayed eigenvalues is remarkably good, so giving ground to the conjecture on the Karcher means reported in the final part of Section 3.2.

4.5. Galerkin Discretization of the Laplacian Eigenvalue Problem

The one-dimensional Laplacian eigenvalue problem is given by the differential equation

$$-u_j^{\prime \prime }\left(x\right)={\lambda }_j{u}_j\left(x\right),x\in (0,1),$$

(22)

with Dirichlet boundary conditions:

$$u_j\left(0\right)=u_j\left(1\right)=0.$$

The goal is to find the eigenvalues ${\lambda }_j\in {\mathbb{R}}^{+}$ and corresponding eigenfunctions $u_j\in H_0^1\left([0,1]\right)$, for $j=1,2,....\infty$, where $H_0^1\left([0,1]\right)$ is the standard Sobolev space of $L^2$ functions with $L^2$ derivatives and vanishing at the boundaries.

Weak formulation

To derive the weak formulation, we multiply both sides of the differential equation by a test function $v\in H_0^1\left([0,1]\right)$ and integrate over the interval $[0,1]$ i.e. $-{\int }_0^1{u}_j^{\prime \prime }\left(x\right)v\left(x\right)dx={\int }_0^1{\lambda }_j{u}_j\left(x\right)v\left(x\right)dx$. Using integration by parts on the left-hand side and noting that the boundary terms vanish, we deduce ${\int }_0^1{u}_j^{\prime \prime }\left(x\right)v\left(x\right)dx={\int }_0^1{u}_j^{\prime }\left(x\right)v^{\prime }\left(x\right)dx$, so that the weak form becomes

$${\int }_0^1{u}_j^{\prime }\left(x\right)v^{\prime }\left(x\right)dx={\lambda }_j{\int }_0^1{u}_j\left(x\right)v\left(x\right)dx,$$

for every test function $v\in H_0^1\left([0,1]\right)$. The latter is rewritten compactly as

$$a(u_j,v)={\lambda }_j(u_j,v),$$

(23)

where the bilinear form $a(u_j,v)$ and the $L^2$ inner product $(u_j,v)$ are defined as

$$a(u_j,v):={\int }_0^1{u}_j^{\prime }\left(x\right)v^{\prime }\left(x\right)dx,(u_j,v):={\int }_0^1{u}_j\left(x\right)v\left(x\right)dx.$$

Galerkin Approximation

The weak formulation allows us to use the Galerkin method to approximate the solution. Let ${\mathcal{W}}_n=\mathrm{span}\{{\varphi }_1,\cdots ,{\varphi }_{N_n}\}$ be a finite-dimensional subspace of $H_0^1\left([0,1]\right)$. The weak problem now becomes: find approximate eigenvalues ${\lambda }_{j,n}\in {\mathbb{R}}^{+}$ and eigenfunctions $u_{j,n}\in {\mathcal{W}}_n$, for $j=1,2....,N_n$ such that, for all $v_n\in {\mathcal{W}}_n$, we have $a(u_{j,n},v_n)={\lambda }_{j,n}(u_{j,n},v_n)$. By expanding $u_{j,n}$ and $v_n$ in terms of the basis functions $\left\{u_i\right\}$, we obtain $u_{j,n}={\sum }_{i=1}^{N_n}{c}_i{\varphi }_i$, $v_n={\sum }_{k=1}^{N_n}{d}_k{\varphi }_k$, and substituting into the bilinear forms, the generalized eigenvalue problem

$$K_n{\mathbf{c}}_j={\lambda }_{j,n}{M}_n{\mathbf{c}}_j$$

is defined, where the stiffness matrix $K_n$ and mass matrix $M_n$ are defined as

$$K_n={\left[a({\varphi }_j,{\varphi }_i)\right]}_{ij=1}^{N_n}=[{\int }_0^1{\varphi }_j^{\prime }\left(x\right){\varphi }_i^{\prime }\left(x\right)dx{]}_{ij=1}^{N_n},$$

(24)

$$M_n={\left[({\varphi }_j,{\varphi }_i)\right]}_{ij=1}^{N_n}=[{\int }_0^1{\varphi }_j\left(x\right){\varphi }_i\left(x\right)dx{]}_{ij=1}^{N_n}.$$

(25)

Both $K_n$ and $M_n$ are symmetric and positive definite, due to the coercive character of the underlying bilinear forms.

4.6. Quadratic C⁰ B-Spline Discretization

In the quadratic C⁰ B-spline discretization of the one-dimensional Laplacian eigenvalue problem, the basis functions ${\varphi }_1,\cdots ,{\varphi }_{N_n}$ are chosen as B-splines of degree 2 defined on a uniform mesh with step size ${\displaystyle \frac{1}{n}}$. The basis functions are explicitly constructed on the knot sequence $\{0,0,0,{\displaystyle \frac{1}{n}},{\displaystyle \frac{1}{n}},{\displaystyle \frac{2}{n}},{\displaystyle \frac{2}{n}}\cdots ,{\displaystyle \frac{n-1}{n}},{\displaystyle \frac{n-1}{n}},1,1,1\},$ (excluding the first and last B-splines, which do not vanish on the boundary of $[0,1]$); see [22]. The resulting normalized stiffness and mass matrices are given by

(26)

(27)

The stiffness matrix ${\displaystyle \frac{1}{n}}{K}_n$ and the mass matrix $n{M}_n$ contain, as principal submatrices, the unilevel 2-block Toeplitz matrices $T_{n-1}\left(f\right)$ and $T_{n-1}\left(h\right)$, according to the notation in Section 2.7 with $d=1$ and $r=2$, where the generating functions $f\left(\theta \right)$ and $h\left(\theta \right)$ are given by

$$f\left(\theta \right):={\displaystyle \frac{1}{3}}\left(\left[\begin{array}{cc}0& -2\\ 0& -2\end{array}\right]e^{i\theta }+\left[\begin{array}{cc}4& -2\\ -2& 8\end{array}\right]+\left[\begin{array}{cc}0& 0\\ -2& -2\end{array}\right]e^{-i\theta }\right),$$

$$f\left(\theta \right)={\displaystyle \frac{1}{3}}\left[\begin{array}{cc}4& -2-2e^{i\theta }\\ -2-2e^{-i\theta }& 8-4cos\theta \end{array}\right],$$

and

$$h\left(\theta \right):={\displaystyle \frac{1}{30}}\left(\left[\begin{array}{cc}0& 3\\ 0& 1\end{array}\right]e^{i\theta }+\left[\begin{array}{cc}4& 3\\ 3& 12\end{array}\right]+\left[\begin{array}{cc}0& 0\\ 3& 1\end{array}\right]e^{-i\theta }\right),$$

$$h\left(\theta \right)={\displaystyle \frac{1}{30}}\left[\begin{array}{cc}4& 3+3e^{i\theta }\\ 3+3e^{-i\theta }& 12+2cos\theta \end{array}\right].$$

Since ${\left\{T_n\left(g\right)\right\}}_n{\sim }_{\mathrm{GLT}}g$ for any Lebesgue integrable generating function g (Axiom GLT 2, part 1), the theory of GLT sequences and a basic use of the extradimensional approach [2,42] lead to

$${\left\{{\displaystyle \frac{1}{n}}{K}_n\right\}}_n{\sim }_{\mathrm{GLT}}f\left(\theta \right),$$

(28)

$${\left\{n{M}_n\right\}}_n{\sim }_{\mathrm{GLT}}h\left(\theta \right).$$

(29)

By the Axiom GLT 3, the linear combination of two GLT sequences is again a GLT sequence, with the symbol being the corresponding linear combination of their symbols. The matrix-sequence ${\left\{L_n\right\}}_n$ with

$$L_n:={\displaystyle \frac{1}{n}}{K}_n+n{M}_n,$$

(30)

is a linear combination of the sequences ${\left\{{\displaystyle \frac{1}{n}}{K}_n\right\}}_n$ and ${\left\{n{M}_n\right\}}_n$. Consequently,

$${\left\{L_n\right\}}_n{\sim }_{\mathrm{GLT}}f\left(\theta \right)+h\left(\theta \right)=e\left(\theta \right),$$

where $f\left(\theta \right)$ and $h\left(\theta \right)$ are the symbols of the stiffness matrix ${\displaystyle \frac{1}{n}}{K}_n$ and the mass matrix $n{M}_n$, respectively.

Figure 5. Comparison between the symbol $G\left(f\right(\theta ),h(\theta ),e(\theta \left)\right)$ and $\mathrm{eig}\left(G(n^{-1}{K}_n,n{M}_n,L_n)\right)$.

Figure 5. Comparison between the symbol $G\left(f\right(\theta ),h(\theta ),e(\theta \left)\right)$ and $\mathrm{eig}\left(G(n^{-1}{K}_n,n{M}_n,L_n)\right)$.

The numerical evidence is convincing: furthermore, from the related figure, we observe that there exist two branches of the spectrum in accordance with Theorem 5 with $r=2$.

4.7. Cubic C¹ B-Spline Discretization

In the cubic C¹ B-spline discretization on a uniform mesh with step size ${\displaystyle \frac{1}{n}}$, the basis functions ${\varphi }_1,\cdots ,{\varphi }_{N_n}$ are chosen as the B-splines; see [22]. The resulting normalized stiffness and mass matrices are given by

(31)

(32)

According to the notation in Section 2.7 with $d=1$ and $r=2$, the stiffness matrix ${\displaystyle \frac{1}{n}}{K}_n$ and the mass matrix $n{M}_n$ contain, as principal submatrices, the unilevel 2-block Toeplitz matrices $T_{n-1}\left(f\right)$ and $T_{n-1}\left(h\right)$, where the generating functions $f\left(\theta \right)$ and $h\left(\theta \right)$ are $2\times 2$ matrix-valued functions given by

$$f\left(\theta \right):={\displaystyle \frac{1}{40}}\left(\left[\begin{array}{cc}-15& -15\\ -3& -15\end{array}\right]e^{i\theta }+\left[\begin{array}{cc}48& 0\\ 0& 48\end{array}\right]+\left[\begin{array}{cc}-15& -3\\ -15& -15\end{array}\right]e^{-i\theta }\right),$$

$$f\left(\theta \right)={\displaystyle \frac{1}{40}}\left[\begin{array}{cc}48-30cos\left(\theta \right)& -15-3e^{-i\theta }\\ -3e^{i\theta }-15e^{-i\theta }& 48-30cos\theta \end{array}\right],$$

and

$$h\left(\theta \right):={\displaystyle \frac{1}{560}}\left(\left[\begin{array}{cc}9& 53\\ 1& 9\end{array}\right]e^{i\theta }+\left[\begin{array}{cc}128& 80\\ 80& 128\end{array}\right]+\left[\begin{array}{cc}9& 1\\ 53& 9\end{array}\right]e^{-i\theta }\right),$$

$$h\left(\theta \right)={\displaystyle \frac{1}{560}}\left[\begin{array}{cc}128+18cos\left(\theta \right)& 80+53e^{i\theta }+e^{-i\theta }\\ 80+53e^{-i\theta }+e^{i\theta }& 128+18cos\left(\theta \right)\end{array}\right].$$

As discussed in Section 4.6, the same reasoning applies to the matrix-sequence ${\left\{L_n\right\}}_n$ with

$$L_n:={\displaystyle \frac{1}{n}}{K}_n+n{M}_n,$$

which results in the GLT sequence symbol $e\left(\theta \right)=f\left(\theta \right)+h\left(\theta \right)$. The numerical evidence is again strong and we can see from the related figure that there exist two branches of the spectrum in accordance with Theorem 5 with $r=2$.

Figure 6. Comparison between the symbol $G\left(f\right(\theta ),h(\theta ),e(\theta \left)\right)$ and $\mathrm{eig}\left(G(n^{-1}{K}_n,n{M}_n,L_n)\right)$.

Figure 6. Comparison between the symbol $G\left(f\right(\theta ),h(\theta ),e(\theta \left)\right)$ and $\mathrm{eig}\left(G(n^{-1}{K}_n,n{M}_n,L_n)\right)$.

4.8. Minimal Eigenvalues and Conditioning

In the last part of the current numerical section, we consider the problem of understanding the extremal behavior of the spectral of the conditioning of ${\left\{G(A_n^{\left(1\right)},A_n^{\left(2\right)},A_n^{\left(3\right)})\right\}}_n$ in dependence of the analytical features of the corresponding GLT symbol. The idea is borrowed from the literature, where we find papers dealing with the extremal eigenvalues in a Toeplitz setting [16,29,30,37], in a r-block Toeplitz setting [31,32], in a differential setting [27,43], including multilevel cases of $d>1$. Here we restrict our attention to the unilevel scalar setting with $d=r=1$ and we consider again the examples in Section 4.1 and in Section 4.3.

4.8.1. Example 1 (1D): Minimal Eigenvalue

We consider the example in Section 4.1. The generating function $\kappa (x,\theta )={(2-2cos\left(\theta \right))}^2$ of $A_n$ has a unique zero of order 4 at $\theta =0$, while the matrix-sequence ${\left\{B_n\right\}}_n$ shows the GLT symbol ${\xi (x,\theta ))=x(2-2cos\left(\theta \right))}^2$ with combinations of zeros at $x=0$ of order 1 and at $\theta =0$ of order 4. According to the theory, the minimal eigenvalue of $A_n$ is positive and tends to zero as $n^{-4}$ (see [16,29,30,37]). On the other hand, according to [27,43], we know that eigenvalue of $B_n$ is positive and tends to zero as $n^{-4}$. If we consider the geometric mean of the two symbols, we deduce that it has again zeros of order at most 4: as a consequence, heuristically we expect that the minimal eigenvalue of $X_n=G(A_n,B_n)$ is positive and tends to zero as $n^{-4}$ as it is perfectly verified in the table below with ${\alpha }_j$ tending to 4.

• $X_n=G(A_n,B_n)$
• Take $n_j=40*2^j,j=0,1,2,3,$
• Compute ${\tau }_j={\lambda }_{min}\left(X_n\right),j=0,1,2,3,$
• Compute ${\alpha }_j={log}_2\left({\displaystyle \frac{{\tau }_j}{{\tau }_{j+1}}}\right),j=0,1,2.$

	${\mathit{\tau }}_{\mathit{j}}$	${\mathit{\alpha }}_{\mathit{j}}$
$n_0=40$	0.00010994	3.8698
$n_1=80$	0.00000752	3.9398
$n_2=160$	0.00000049	4.0297
$n_3=320$	0.00000003

	${\mathit{\tau }}_{\mathit{j}}$	${\mathit{\alpha }}_{\mathit{j}}$
$n_0=40$	0.00010994	3.8698
$n_1=80$	0.00000752	3.9398
$n_2=160$	0.00000049	4.0297
$n_3=320$	0.00000003

4.8.2. Example 3 (1D): Minimal Eigenvalue

We consider the example in Section 4.3. The generating function ${\kappa }_1$ of $A_n^{\left(1\right)}$ is strictly positive, the generating function ${\kappa }_2(x,\theta )=x^2$ of $A_n^{\left(2\right)}$ has a unique zero of order 2 at $x=0$, while the matrix-sequence ${\left\{A_n^{\left(3\right)}\right\}}_n$ shows the GLT symbol ${\kappa }_3(x,\theta ))=x^4{(2-2cos\left(\theta \right))}^4$ with combinations of zeros at $x=0$ of order 4 and at $\theta =0$ of order 4. According to the theory, the minimal eigenvalue of $A_n^{\left(1\right)}$ is positive and tends to $1=min{\kappa }_1$ as (see [16,29,30,37]). On the other hand, according to [27,43], we know that eigenvalue of $A_n^{\left(3\right)}$ is positive and tends to zero as $n^{-4}$, while for $A_n^{\left(2\right)}$ it is trivial tp check that the minimal eigenvalue tends to zero as $n^{-2}$. If we consider the geometric mean of the two symbols, we observe that it has zeros of order at most 2: hence, heuristically we expect that the minimal eigenvalue of $X_n=G(A_n^{\left(1\right)},A_n^{\left(2\right)},A_n^{\left(3\right)})$ is positive and tends to zero as $n^{-2}$ as it is perfectly verified in the table below with ${\alpha }_j$ tending to 2.

• $X_n=G(A_n^{\left(1\right)},A_n^{\left(2\right)},A_n^{\left(3\right)})$
• Take $n_j=40*2^j,j=0,1,2,3,$
• Compute ${\tau }_j={\lambda }_{min}\left(X_n\right),j=0,1,2,3,$
• Compute ${\alpha }_j={log}_2\left({\displaystyle \frac{{\tau }_j}{{\tau }_{j+1}}}\right),j=0,1,2.$

	${\mathit{\tau }}_{\mathit{j}}$	${\mathit{\alpha }}_{\mathit{j}}$
$n_0=40$	0.0014	2.2223
$n_1=80$	0.00034797	2
$n_2=160$	0.000086993	2
$n_3=320$	0.000021748

	${\mathit{\tau }}_{\mathit{j}}$	${\mathit{\alpha }}_{\mathit{j}}$
$n_0=40$	0.0014	2.2223
$n_1=80$	0.00034797	2
$n_2=160$	0.000086993	2
$n_3=320$	0.000021748

5. Conclusions

We have studied the spectral distribution of the geometric mean of two or more matrix-sequences constituted by HPD matrices, under the assumption that k input matrix-sequences belong to the same d-level, r-block GLT *-algebra, where $k\ge 2$, $d,r\ge 1$. For $k=2$ an explicit formula exists and this has allowed to prove that the new matrix-sequence is of GLT type for any fixed $d,r$, with GLT symbol being the corresponding geometric mean of the $k=2$ input symbols, according to Theorems 4 and 5. As a consequence the spectral distribution of the geometric mean of $k=2$ matrix-sequences is formally demonstrated. For $k>2$, an explicit formula is not available and we have considered the Karcher mean, for which efficient iterative procedures available in the literature can be used for computational purposes. Interestingly enough, using the same tools as in the proof of Theorems 4 and 5, we have shown that all the matrix-sequences made by the Karcher mean iterates are still GLT matrix-sequences with symbols converging the geometric mean of the k symbols, provided that the matrix-sequence of the initial input matrices is a GLT matrix-sequence made by HPD matrices. The theoretical results have been confirmed through various numerical experiments, where we compared the eigenvalues of the geometric mean with a uniform sampling of the geometric mean of the symbols. A very good agreement has been observed in all the numerical tests and aven for very moderate matrix-sizes. In fact, as the matrix-size increases, the symbol provides a better and better approximation of the eigenvalues, both in one-dimensional and two-dimensional cases, and both with scalar and block matrices, i.e. using GLT matrix-sequences taken from the differential world with $k=2,3$, $d,r=1,2$ and including finite difference, finite element, isogeometric analysis approximations of constant and variable differential operators. Few open questions remain and in the subsequent lines we indicate three main issues:

• A formal proof of the GLT nature of the Karcher mean of $k>2$ HPD GLT matrix-sequences has to be given under the assumption that the initial guess is a HPD GLT matrix-sequence; in this respect, also from a computational viewpoint, starting with a initial guess ${\left\{X_{\mathbf{n},0}\right\}}_{\mathbf{n}}$ having already as GLT symbol the geometric mean of the input symbols should reduce sensibly the number of iterations;
• in connection with Theorems 4 and 5 and $k=2$, the ALM axioms suggest that the technical assumption on the invertibility almost everywhere of the input GLT symbols is not necessary for every $d,r\ge 1$;
• a completely open problem concerns the study of the extremal eigenvalues of the geometric means of GLT matrix-sequences as a function of the analytic features of the geometric means of the GLT symbols: in this direction it should be recalled that a rich literature exists regarding the extremal eigenvalues in a Toeplitz setting [16,29,30,37], in a r-block Toeplitz setting [31,32], in a differential setting [27,43], so involving all the types of examples considered in the numerical experiments. Prelinary numerical experiments in the unilevel scalar setting with $k=2,3$ has been performed in Section 4.8 and the results are quite promising. Interestingly enough, and substantially mimicking the cases already studied in the literature, it seems that the order of the zeros of the GLT symbol decides the asymptotic behavior of the minimal eigenvalues and hence of the conditioning of $G(A_n^1,\dots ,A_n^k))$, at least for $d=r=1$, $k=2,3$.

All the previous items show that a lot of scientific research has still to be done in connection with the findings discussed in the current work.