Krylov Iterative Methods for the Geometric Mean of Two Matrices Times a Vector

TL;DR

This paper introduces an efficient Krylov space-based method for computing the geometric mean of two positive definite matrices times a vector, leveraging matrix-vector products for low-cost approximations.

Contribution

It presents a novel Krylov iterative approach specifically designed for the geometric mean of matrices times a vector, improving computational efficiency.

Findings

Method achieves accurate approximations with fewer matrix-vector products

Applicable to large-scale positive definite matrices

Reduces computational cost compared to traditional methods

Abstract

In this work, we are presenting an efficient way to compute the geometric mean of two positive definite matrices times a vector. For this purpose, we are inspecting the application of methods based on Krylov spaces to compute the square root of a matrix. These methods, using only matrix-vector products, are capable of producing a good approximation of the result with a small computational cost.

Tables

Table 1. Table 1: Time comparisons between the proposed methods.

Method	Measured time
Dimension	$1600\times 1600$	$2500\times 2500$	$3600\times 3600$	$4900\times 4900$
Minimax	42.3826	136.398	382.048	743.067
Schur	24.3768	90.8797s	253.466	651.758
Leja Points	2.26687	3.71883s	6.18135	9.54854
Adaptive	2.03123	3.54691	5.5625	9.01572

Equations

$$\underset{n,m\to \infty }{\lim }\frac{\mathcal{S}(n,m)}{nm}=0.$$

$$A(A^{-1}B{)}^{1/2}.$$

$$M^{1/2}=K\text{mbox}diag(\sqrt{{\lambda }_1},\dots ,\sqrt{{\lambda }_n})K^{-1}$$

$$A\#B=A(A^{-1}B{)}^{1/2}=(B{A}^{-1}{)}^{1/2}A=B(B^{-1}A{)}^{1/2}B=B\#A.$$

$$f(z)={\int }_{\Gamma }\frac{d\gamma (x)}{z-x},$$

$$f(z)=z^{-1/2}={\int }_{-\infty }^0\frac{1}{z-x}\frac{dx}{\pi \sqrt{-x}}.$$

$$K_m(A,b)=\text{mbox}span\{b,Ab,A^{2}b,\dots ,A^{m-1}b\}\subseteq {\mathbb{C}}^n.$$

$$\mbox{dim }K_{i}(A,b)=\left\{\begin{array}[]{ll}i&\mbox{if }i<k,\\ k&\mbox{otherwise}.\end{array}\right.$$

$$A{Q}_k=B{Q}_k{T}_k+{\beta }_{k+1}B{q}_{k+1}{e}_k^T{Q}_k^{\ast }B{Q}_k=I_k,$$

$$Q^{\ast }AQ=T,Q^{\ast }BQ=I.$$

$$B\#A=B(B^{-1}A{)}^{1/2}=(Q{T}^{-1/2}{Q}^{\ast }{)}^{-1}=BQ{T}^{1/2}{Q}^{\ast }B.$$

$$(B\#A)v\approx B{Q}_k{T}_k^{1/2}{e}_1\Vert v{\Vert }_B,$$

$$Q^{\ast }AQ=H,$$

$$A{q}_k=\sum _{i=1}^{k+1}{h}_{ik}{q}_i,\text{mbox}withk=1,\dots ,n-1.$$

$$h_{k+1,k}{q}_{k+1}=A{q}_k-\sum _{i=1}^k{h}_{ik}{q}_i=r_k,$$

$$h_{ik}=q_i^{\ast }A{q}_k\text{mbox}withi=1,\dots ,k.$$

$$\text{mbox}span\{q_1,q_2,\dots ,q_k\}=\text{mbox}span\{q_1,A{q}_1,\dots ,A^{k-1}{q}_1\},$$

$$A{Q}_k=Q_k{H}_k+h_{k+1,k}{q}_{k+1}{e}_k^T,$$

$$Q_k^{\ast }A{Q}_k=H_k,$$

$$f_k=Q_{k}f(H_k)e_1\Vert b{\Vert }_2=Q_{k}f(H_k)Q_k^{\ast }b,$$

$$(A\#B)v\approx A{Q}_k{H}_k^{-1/2}{Q}_k^{\ast }v=A{Q}_k{H}_k^{-1/2}{e}_1\Vert v{\Vert }_2.$$

$$q_{m-1}(z)=\prod _{j=1}^{m-1}(1-z/{\xi }_j),$$

$$Q_m(A,b)=q_{m-1}(A{)}^{-1}\text{mbox}span\{b,Ab,\dots ,A^{m-1}b\}\subseteq {\mathbb{C}}^n.$$

$$x_j=(I-A/{\xi }_j{)}^{-1}A{v}_j$$

$$x_j=\sum _{i=1}^{j+1}{v}_i{h}_{i,j},$$

$$A\left(v_j+\sum _{i=1}^{j+1}{v}_i{h}_{i,j}{\xi }_j^{-1}\right)=\sum _{i=1}^{j+1}{v}_i{h}_{i,j},$$

$$A{V}_m(I_m+H_m{D}_m)+A{v}_{m+1}{h}_{m+1,m}{\xi }_m^{-1}{e}_m^T=V_m{H}_m+v_{m+1}{\xi }_m^{-1}{e}_m^T,$$

$$\underline{H_{m}}=\left[\begin{array}[]{c}H_{m}\\ h_{m+1,m}e_{m}^{T}\end{array}\right]\qquad\mbox{and}\qquad\underline{K_{m}}=\left[\begin{array}[]{c}I_{m}+H_{m}D_{m}\\ h_{m+1,m}{\xi}_{m}^{-1}e_{m}^{T}\end{array}\right]$$

$$A{V}_{m+1}\underline{K_m}=V_{m+1}\underline{H_m},$$

$$A{V}_m{K}_m=V_{m+1}\underline{H_m},$$

Full text

K****eywords geometric mean  $\cdot$ Krylov spaces  $\cdot$ Krylov methods  $\cdot$ iterative methods  $\cdot$ positive definite matrices  $\cdot$ sparse matrices  $\cdot$ rational Arnoldi method

1 Introduction

Many problems in engineering, medicine and computer science make use of the geometric mean of two matrices $A\#B=A(A^{-1}B)^{1/2}$, with $A$ and $B$ being positive definite matrices. These applications include the calculation of electrical networks [1], diffusion tensor imaging [17] and image deblurring [9]. For these applications there are methods to approximate a function of a matrix, as shown in [15] and [13], that are often used with acceptable results in term of computational cost.

But in the solution of elliptic partial differential equations with domain decomposition (such as in [2] and [3]), in certain cases it is required to compute the geometric mean of two matrices times a vector $(A\#B)v$, where $A$ and $B$ are often very large and sparse. A matrix of dimension $n\times m$ is called sparse if defined $\mathcal{S}(n,m)$ as the number of its non-zero elements, it holds

[TABLE]

The use of algorithms created to compute $A\#B$ could lead to unsustainable computational times or sometimes lead to the impossibility of calculating the result, since these methods do not exploit the sparse structure of data, that instead should represent an advantage and a saving in terms of operations.

For these reasons, methods based on Krylov spaces are preferred. In [2] and [18], the usage of a generalized version of the Lanczos method is exploited (we briefly present it in Section 3). In this work, we will use methods based on both polynomial and rational Krylov spaces in order to compute the geometric mean of two matrices times a vector in an efficient way, by exploiting the sparse structure of the involved matrices. In fact, these methods are easier to use, because they only compute matrix-vector products and solve linear systems with the conjugate gradient method. Moreover, their use as iterative methods is another advantage to reduce the computational cost. Specifically, we will focus on the Arnoldi method and on some variations of the rational Arnoldi method, explained in Section 3 and Section 4 respectively. Experiments presented in Section 5 show how these methods, especially rational Arnoldi ones, are efficient in terms of approximation and computational time.

2 Geometric Mean of Two Matrices

The geometric mean $A\#B$ of two positive definite matrices $A,B\in\mathbb{C}^{n\times n}$ is defined in [5] as

[TABLE]

Let $M$ be a diagonalizable matrix with positive eigenvalues, so that an invertible matrix $K$ and a diagonal matrix $D=\mbox{diag}({\lambda}_{1},\ldots,{\lambda}_{n})$ exist such that $M=KDK^{-1}$. The square root of M is defined as

[TABLE]

that turns out to be still positive definite. The matrix $A^{-1}B$ is diagonalizable with positive eigenvalues, because it is similar to $(A^{1/2})^{-1}B(A^{1/2})^{-1}$ which is positive definite, so the notation $A(A^{-1}B)^{1/2}$ makes sense. It can be shown that $(A^{1/2})^{-1}=(A^{-1})^{1/2}$, then we use the notation $A^{-1/2}$ without ambiguity. We observe also that

[TABLE]

We also have that, if both $A$ and $B$ are positive definite matrices, then so is $A(A^{-1}B)^{1/2}$.

To get an efficient computation of $(A\#B)v$ it is important not to use expensive operations. The most problematic part is to find a good approximation of $A(A^{-1}B)^{1/2}v=A(B^{-1}A)^{-1/2}v$ without explicitly computing $B^{-1}$ or the square root of $B^{-1}A$. We recall the fact that $z^{-1/2}$ is a Markov function, i.e. a function of the form

[TABLE]

with $\gamma$ complex measure over the closed set $\Gamma\subseteq\mathbb{C}$, because

[TABLE]

3 Polynomial Krylov Spaces

Given $A\in\mathbb{C}^{n\times n}$ and $b\in\mathbb{C}^{n}$, the $m$-th Krylov polynomial space $K_{m}(A,b)$ associated with them is defined as

[TABLE]

With the increase of $m$, we get Krylov spaces all nested one inside another, $K_{1}(A,b)\subseteq K_{2}(A,b)\subseteq\ldots\subseteq K_{d}(A,b)$. It can be shown that if $b\neq 0$, there exists $k\leq n$ such that

[TABLE]

For each step the space $K_{i}(A,b)$ has a basis $\{b,Ab,\ldots,A^{i-1}b\}$, but from a certain $k+1$ on the new term $A^{k}b$ is linearly dependent with the other vectors $b,\ldots,A^{k-1}b$ in the basis, giving that $K_{k+1}(A,b)$ is the same of $K_{k}(A,b)$.

3.1 Generalized Lanczos Method

First we briefly recall the variation of the Lanczos method presented in [18] and [2], called the generalized Lanczos method (see [18], [6] and [10] for an introduction to the Lanczos method). This iterative method is adapted in [2] in order to compute the geometric mean of two positive definite matrices times a vector. Given two positive definite matrices $A,B\in\mathbb{C}^{n\times n}$, the method constructs at the step $k$ a new vector $q_{k}\in\mathbb{C}^{n}$ that is $B$-orthogonal, that is $q_{i}^{*}Bq_{j}=0$ with $i\neq j$, such that for each step we have

[TABLE]

with $Q_{k}=[q_{1},\ldots,q_{k}]$, $I_{k}$ the $k\times k$ identity matrix and $e_{k}$ its $k$-th column.

The matrix $T_{k}\in\mathbb{C}^{k\times k}$ is the projection of $A$ into the space generated by the columns of the matrix $Q_{k}$. At the $n$-th step this gives us the equalities

[TABLE]

Using this method we find an expression for $B\#A$ (that is the same of $A\#B$) as

[TABLE]

To compute $(B\#A)v$, with $v\in\mathbb{R}^{n}$, we can stop at the $k$-th step, and taking $q_{1}=v/\|v\|_{B}$, we have

[TABLE]

because $Q_{k}^{*}Bv=e_{1}\|v\|_{B}$, that is a better approximation the closer $k$ is to $n$.

3.2 Arnoldi Method

The Arnoldi method, exposed in [14], is a general version of the Lanczos method that does not need the starting matrices to be positive definite. This method computes the decomposition

[TABLE]

where $Q\in\mathbb{C}^{n\times n}$ is unitary, and the vectors $q_{1},\ldots,q_{n}$ are its columns, and $H\in\mathbb{C}^{n\times n}$ is upper Hessenberg. From the equality (13) we get that

[TABLE]

The previous relation can be rewritten as

[TABLE]

and, because $Q$ is unitary, we have that

[TABLE]

If $r_{k}\neq 0$, so $q_{k+1}=r_{k}/h_{k+1,k}$, where $h_{k+1,k}=\|r_{k}\|_{2}$. The Arnoldi method is a Krylov space method because

[TABLE]

that is the vectors $q_{1},\ldots,q_{k}$ constructed by the Arnoldi method are a basis for the Krylov space $K_{k}(A,b)$. Used as an iterative method, at the $k$-th step the Arnoldi method give us the factorization

[TABLE]

where $Q_{k}$ consist of the first $k$ column of $Q$, $H_{k}$ is a $k\times k$ upper Hessenberg matrix and $e_{k}$ is the $k$-th column of the $k\times k$ identity matrix. Because

[TABLE]

we have that $H_{k}$ is the projection of $A$ in $K_{k}(A,b)$. The orthonormalization of vectors $q_{1},\ldots,q_{n}$ comes from the modified Gram-Schmidt method (see [22]), but in finite precision arithmetic a loss of orthogonality can occur, as shown in [18].

Using the Arnoldi method we can approximate $f(A)b$, starting with the vector $q_{1}=b/\|b\|_{2}$. As shown in [14] and [4] it stands that, at the $k$-th step, the approximation is

[TABLE]

that is the same of calculating $f$ into the Krylov space $K_{k}(A,B)$ and then expanding the result into the original space $\mathbb{C}^{n\times n}$. So an approximation for $(A\#B)v=A(A^{-1}B)^{1/2}v$ is

[TABLE]

4 Rational Krylov Spaces

The definition of rational Krylov spaces is similar to the polynomial one, except for the presence of a denominator (see [11] and [12] for an introduction). Given $A\in\mathbb{C}^{n\times n}$ and $b\in\mathbb{C}^{n}$, and the sequence of polynomials $q_{m-1}(A)$ of degree $m-1$

[TABLE]

in which the values ${\xi}_{j}\in\mathbb{C}\cup\{\infty\}$ are called poles and are numbers in the extended complex plane different from all the eigenvalues of $A$ and [math], the rational Krylov space of order $m$ associated to them is defined as

[TABLE]

In the case in which ${\xi}_{j}=\infty$ for some $j$ the corresponding factor $(1-z/{\xi}_{j})$ in (23) is replaced by $1$.

In the previous definition it is not taken into account the case ${\xi}_{j}=0$, but we can exclude every other value $\sigma$ simply by using the new values $\widehat{A}=A-\sigma I$ and $\widehat{{\xi}_{j}}={\xi}_{j}-\sigma$. This can be done considering the new polynomial $q_{m-1}$ that uses poles $\widehat{\xi}_{j}$ and the new matrix $\widehat{A}$ instead of $A$. Like the polynomial ones, rational Krylov spaces are also of increasing dimension and are nested until a certain dimension $k$ from which they do not change anymore. Polynomial Krylov spaces can be seen as a special case of the rational ones, when every pole ${\xi}_{j}$ is $\infty$, that is when $q_{m-1}=1$.

4.1 Rational Arnoldi Method

To compute a basis $V_{m}=[v_{1},\ldots,v_{m}]\in\mathbb{C}^{n\times m}$ for the rational Krylov space $Q_{m}(A,b)$ we can use the rational Arnoldi method, as proposed in [11] and in [4]. Starting with $v_{1}=b/\|b\|_{2}$, in the following iterations the vector $v_{j+1}$ is generated orthogonalizing

[TABLE]

against the previous orthonormal vectors $v_{1},\ldots,v_{j}$. If $h_{j+1,j}\neq 0$, we have that

[TABLE]

and $v_{j+1}=x_{j}/h_{j+1,j}$, instead if $h_{j+1,j}=0$ we can proceed finding a new vector $v_{j+1}$ that is orthonormal to $v_{1},\ldots,v_{n}$.

We also have that

[TABLE]

that gives us the decomposition

[TABLE]

in which $D_{m}=\mbox{diag}({\xi}_{1}^{-1},\ldots,{\xi}_{m}^{-1})$, $I_{m}$ is the $m\times m$ identity matrix and $e_{m}$ is its $m$-th column. Defining

[TABLE]

we have that

[TABLE]

where $V_{m+1}=[V_{m},v_{m+1}]$. Finally, if the last pole ${\xi}_{m}$ is $\infty$, the decomposition is simplified into

[TABLE]

where $K_{m}$ is the $m\times m$ upper part of $\underline{K_{m}}$.

Computed the basis $V_{m}$ of $Q_{m}(A,b)$, we can approximate $f(A)b$ as

[TABLE]

The benefit of this method is that, in many interesting cases, $f_{m}^{RA}$ is a good approximation of $f(A)b$ also if $m$ is small: it is required to compute $f(A_{m})$, but $A_{m}$ is small compared to $A$ itself. If the last pole ${\xi}_{m}$ is $\infty$, it is not even necessary to compute $A_{m}=V_{m}^{*}AV_{m}$, but it can be computed $A_{m}=H_{m}K_{m}^{-1}$. Also $H_{m}$ and $K_{m}$ are small compared to $A$, so computing $K_{m}^{-1}$ or the product $H_{m}K_{m}^{-1}$ is possible. For our problem $(A\#B)v$ we have found the approximation

[TABLE]

starting the method with $v_{1}=v/\|v\|_{2}$.

From a Crouzeix theorem exposed in [7] we know that there exists a constant $C\leq 11.08$ for which $\|f(A)\|_{2}\leq C\|f\|_{\Sigma}$, where the second norm indicates the maximum absolute value of $f$ over the compact set $\Sigma$ on which the function is approximated, that in this case is a set containing the spectrum of $A$. As explained in [11], this helps us to find a result of semi-optimality for the approximation given by the rational Arnoldi method: in fact, if $f$ is analytic in $\Sigma$ and we define $f_{m}=r_{m}(A)b$ with $r_{m}\in\mathcal{P}_{m-1}/q_{m-1}$, where $\mathcal{P}_{m-1}$ is a polynomial of degree $m-1$, so it stands

[TABLE]

It seems obvious that the choice of the poles ${\xi}_{j}$ is strongly connected to the function $f$ which we have to approximate. To find a small value of $\|f(A)b-f_{m}^{RA}\|_{2}$, it is necessary to search a sufficiently uniform approximation of $f$ over $\Sigma$. Moreover, sometimes it is convenient to consider only a certain set $\Xi$ for the poles disjoint from $\Sigma$. For Markov functions, and especially for the $z^{-1/2}$ function that appears in our problem, a good set $\Xi$ is $[-\infty,0]$, and as $\Sigma$ we can use $[{\lambda}_{min},{\lambda}_{max}]$, with ${\lambda}_{min}$ and ${\lambda}_{max}$ respectively the minimum and maximum eigenvalues of $A$ (as above, see [11] for further informations). Computing the eigenvalues of $A$ can be computationally expensive, however only the extreme eigenvalues are required.

4.1.1 Extended Krylov Method

A special case of the rational Arnoldi method is the so-called extended Krylov method, proposed in [8] and then investigated in [21] and [16]. It simply chooses alternately the poles as ${\xi}_{even}=\infty$ and ${\xi}_{odd}=0$. This method is good for approximate Markov functions. Its benefit is obvious: the choice of the poles is completely a priori, without information on the spectrum of $A$. It had been observed in [8] that in some cases this is equivalent to a rational Arnoldi method in which all the poles are chosen as an unique asymptotically optimal value. However, finding the optimal pole needs some informations on the spectrum of $A$.

4.1.2 Generalized Leja Points

An alternative way to choose poles, exposed in [11], is the so-called generalized Leja points method, that rely on some logarithmic potential theory’s instruments (for a complete mathematical explanation about it see [19] and [20]). Given two closed sets $\Sigma$ and $\Xi$, both of non-zero logarithmic capacity (the logarithmic capacity of a set is defined as $\mbox{cap}(E)=e^{V_{E}}$, where $V_{E}$ is the Robin constant for $E$) and of positive distance, the pair $(\Sigma,\Xi)$ is called a condenser, and we can associate a number to it called condenser capacity $\mbox{cap}(\Sigma,\Xi)$ (still see [20] for further explanations). We now consider the sequence of functions

[TABLE]

with nodes ${\sigma}_{j}\in\Sigma$ and poles ${\xi}_{j}\in\Xi$. Our aim is to make this sequence in absolute value as large as possible over $\Xi$ and as small as possible over $\Sigma$. It can be shown that for this kind of sequences the relation

[TABLE]

lasts.

To find sequences of functions in which this inequality is true is called the generalized Zolotarev problem for the condenser $(\Sigma,\Xi)$.

A practical method to obtain such functions is the following greedy algorithm: starting with ${\sigma}_{1}$ and ${\xi}_{1}$ as points of minimum distance from $\Sigma$ and $\Xi$, the following points ${\sigma}_{j+1}$ and ${\xi}_{j+1}$ are recursively determined in a way such that lasts

[TABLE]

The points $\{({\sigma}_{j},{\xi}_{j})\}$ are called generalized Leja points.

We want to approximate the Markov function $z^{-1/2}$. In order to archive this, in [11] it is suggested to choose sets $\Sigma$ and $\Xi$ as $\Sigma=[{\lambda}_{min},{\lambda}_{max}]$ and $\Xi=[-\infty,0]$. These sets are optimal respectively for the choice of nodes and poles for that function.

4.1.3 Adaptive Poles

To approximate Markov functions, another way to choose the poles for the rational Arnoldi method is also qsuggested in [12]: they are called adaptive poles. We consider the function

[TABLE]

in which ${\sigma}_{k}$ are the Ritz values of the projection of $A$ at the $m$-th step, i.e. the eigenvalues of $A_{m}$, and $q_{m-1}$ is the denominator associated with the corresponding rational Krylov space. Our aim is to minimize the error in the approximation for each step, and to do so we have to make $|s_{m}(z)|$ uniformly large on the set $\Xi$ choosing the next pole ${\xi}_{m}$ as

[TABLE]

This method is substantially a black-box method, because we do not need any information on the spectrum of $A$, but we only have to compute the eigenvalues of the matrix $A_{m}$, that are considerably less than the eigenvalues of $A$ itself.

5 Experiments and Conclusions

We are now going to consider some experiments in order to test the accuracy of the methods presented above. The first important thing is to analyse their convergence rate, that is how close they are to approximate the value of $(A\#B)v$. For this experiment we have used two positive definite matrices $A\in\mathbb{R}^{100\times 100}$ and $B\in\mathbb{R}^{100\times 100}$ and a vector $v\in\mathbb{R}^{100}$ randomly generated. We have run the methods up to 30 steps for each one.

The results are shown in the following graph, where we have on the x-axis the number of steps (or the number of poles considered for the Minimax method and the Gauss-Chebyshev quadrature, explained in [22] and [15], respectively) and on the y-axis the relative errors with respect to the real value of $(A\#B)v$ in a logarithmic scale.

Figure 1 shows how some methods, like the rational Arnoldi method with Leja points or adaptive poles, behave similarly to the Minimax method (they have almost the same convergence rate, with peaks in their graphics due to some numerical instabilities), but for the latter the knowledge of the full spectrum of the matrices is required.

Another important aspect we had to consider is how much time is required by these methods to compute the results. In fact, taking advantage of the sparse structure, Krylov methods use less operations than those methods that do not exploit it. We have used matrices $A$ and $B$ structured respectively as the finite-difference discretization matrices of the 1D and 2D Laplacian of variable dimension and a vector $v$ with every component equals to $1$. Both matrices are implemented using the Octave CSR sparse format.

[TABLE]

We have compared various methods: a method based on Schur decomposition presented in [15] (that directly compute $A\#B$ and only later makes the product with $v$), the Minimax method presented in [13] (based on contour integrals) and two versions of the rational Arnoldi method, the one with Leja points and the adaptive one. Every method is run for 30 steps. We have measured the execution time of these methods with the Octave commands tic and toc (we have considered the time required to compute the extreme eigenvalues of $B^{-1}A$ with the Octave function eigs() in the estimation for the rational Arnoldi method with Leja points).

Table 1 shows how the computational time required by the Minimax method and the method based on Schur decomposition grows very quickly. For matrices of dimension $10000\times 10000$, the estimated times for the aforementioned methods become practically unsustainable, while methods based on Krylov spaces require only few seconds. These results show how important it is to preserve and exploit the sparse structure of the starting matrices in order to complete these computations in a reasonable time, gaining in terms of operations.