Geometrical Inverse Preconditioning for Symmetric Positive Definite Matrices

For the numerical solution of (3), we start by considering the classical gradient iterations that, from an initial guess $X_0$, are given by where ${\alpha }_k>0$ is a suitable step length. A standard approach is to use the optimal choice i.e., the positive step length that (exactly) minimizes the function $F(X)$ along the negative gradient direction. We present a closed formula for the optimal choice of step length in a more general setting, assuming that the iterative method is given by: where $D_k$ is a search direction in the space of matrices.

Differentiating $\psi (\alpha )$, using that $\langle X^{(k)}A,X^{(k)}A\rangle =n$, and also that and then forcing ${\psi }^{\prime }(\alpha )=0$ the result is obtained, after some algebraic manipulations. ☐

Since ${\parallel I\parallel }_F=\sqrt{n}$, the gradient iterations can be written as which can be further simplified by imposing the condition for uniqueness $\parallel X^{(k)}{A\parallel }_F=\sqrt{n}$. In that case, we set and then we multiply the matrix $Z^{(k+1)}$ by the factor $\sqrt{n}/{\parallel Z^{(k+1)}A\parallel }_F$ to guarantee that $X^{(k+1)}\in S$, i.e., such that $\parallel X^{(k+1)}{A\parallel }_F=\sqrt{n}$.

Concerning the condition that the sequence $\{X^{(k)}\}$ remains in T, in our next result, we establish that if the step length ${\alpha }_k$ remains uniformly bounded from above, then $trace(X^{(k)}A)>0$ for all k.

It follows that and so

Now, since $trace(A^2)=\langle A,A\rangle ={\parallel A\parallel }_F^2$ and we obtain that

Since $0<{\alpha }_k\le \frac{n^{3/2}}{{\parallel A\parallel }_F^2}$, then $(1-\frac{{\alpha }_k}{n^2}\sqrt{n}\parallel A{\parallel }_F^2)>0$, and we conclude that

Since $X^{(k+1)}$ is obtained as a positive scaling factor of $Z^{(k+1)}$, then $w_{k+1}>0$, and the result is established. ☐

Now, for some given matrices A, we cannot guarantee that the step length computed as the absolute value of (4) will satisfy ${\alpha }_k\le (n^{3/2})/{\parallel A\parallel }_F^2$ for all k. Therefore, if $trace(X^{(k+1)}A)=\langle X^{(k+1)}A,I\rangle <0$, then we will set in our algorithm $X^{(k+1)}=-X^{(k+1)}$ to guarantee that $trace(X^{(k+1)}A)\ge 0$, and hence that the cosine between $X^{(k+1)}A$ and I is nonnegative, which is a necessary condition to guarantee that $X^{(k+1)}$ remains in the $PSD$ cone (see, e.g., [24,26,27].

We now present our steepest descent gradient algorithm that will be referred as the CauchyCos Algorithm.

We note that if we start from $X^{(0)}$ such that $\parallel X^{(0)}{A\parallel }_F=\sqrt{n}$ then by construction $\parallel X^{(k)}{A\parallel }_F=\sqrt{n}$, for all $k\ge 0$; for example, $X^{(0)}=(\sqrt{n}/{\parallel A\parallel }_F)I$ is a convenient choice. For that initial guess, $trace(X^{(0)}A)=\langle X^{(0)}A,I\rangle >0$ and again by construction all the iterates will remain in the $PSD$ cone. Notice also that, at each iteration, we need to compute the three matrix–matrix products: $X^{(k)}A$, $(\frac{w_k}{n}{X}^{(k)}A-I)A$, and $\nabla F(X^{(k)})A$, which, for dense matrices, require $n^3$ floating point operations (flops) each. Every one of the remaining calculations (inner products and Frobenius norms) are obtained with n column-oriented inner products that require n flops each. Summing up, in the dense case, the computational cost of each iteration of the CauchyCos Algorithm is $3n^3+O(n^2)$ flops. In Section 2.5, we will discuss a sparse version of the CauchyCos Algorithm and its computational cost.

In spite of the resemblance with the classical Cauchy method for convex constrained optimization, the CauchyCos algorithm involves certain key ingredients in its formulation that splits it apart from the Cauchy method. Therefore, a specialized convergence analysis is required. In particular, we note that the constraint set S on which the iterates are computed, thanks to the scaling step 7, is not a convex set.

We start by establishing the commutativity of all iterates with the matrix A.

It is worth noticing that, using Lemma 2.6 and (5), it follows by simple calculations that $Z^{(k)}$ as well as $X^{(k)}$ are symmetric matrices for all k. In turn, if $X^{(0)}A=A{X}^{(0)}$, this clearly implies using Lemma 2.6 that $X^{(k)}A$ is also a symmetric matrix for all k.

Our next result establishes that the sequences generated by the CauchyCos Algorithm are uniformly bounded away from zero, and hence the algorithm is well-defined.

In that case, the negative gradient, $-\nabla F(\widehat{X})\ne 0$, is a descent direction for the function F at $\widehat{X}$. Hence, there exists $\widehat{\alpha }>0$ such that

Consider now an auxiliary function $\theta :R^{n\times n}\to R$ given by

Clearly, θ is a continuous function, and then $\theta (X^{(k_j)})$ converges to $\theta (\widehat{X})=\delta$. Therefore, for all $k_j$ sufficiently large,

Now, since ${\alpha }_{k_j}$ was obtained using Lemma 2.4 as the exact optimal step length along the negative gradient direction, then using Remark 2.1, it follows that and thus, for all $k_j$ sufficiently large.

On the other hand, since F is continuous, $F(X^{(k_j)})$ converges to $F(\widehat{X})$. However, the whole sequence $\left\{F(X^{(k)})\right\}$ generated by the CauchyCos Algorithm is decreasing, and so $F(X^{(k)})$ converges to $F(\widehat{X})$, and since F is bounded below then for $k_j$ large enough which contradicts (6). Consequently, $\nabla F(\widehat{X})=0$.

Now, using Lemma 2.1, it follows that $\nabla F(\widehat{X})=0$ implies $\widehat{X}=A^{-1}$. Hence, the subsequence $\left\{X^{(k_j)}\right\}$ converges to $A^{-1}$. Nevertheless, as we argued before, the whole sequence $F(X^{(k)})$ converges to $F(A^{-1})=0$, and by continuity the whole sequence $\left\{X^{(k)}\right\}$ converges to $A^{-1}$. ☐

Notice that Theorem 2.2 states that the sequence $X^{(k)}$ converges to $A^{-1}$ which is in the $PSD$ cone. Hence, if $X^{(0)}$ is symmetric and $X^{(0)}A=A{X}^{(0)}$, then $X^{(k)}A$ is symmetric (Lemma 2.6), and after a k iterations (k large enough), the eigenvalues of $X^{(k)}A$ are strictly positive: as a consequence, $X^{(k)}A$ is in the $PSD$ cone.

To avoid the zig-zagging trajectory of the optimal gradient iterates, we now consider a different search direction: to move from $X^{(k)}\in S\cap T$ to the next iterate. Notice that ${\widehat{D}}_{k}A=-\nabla F(X^{(k)})$ and so ${\widehat{D}}_k$ can be viewed as a simplified version of the search direction used in the classical steepest descent method. Moreover, the direction ${\widehat{D}}_k$ in (7) is obtained from the direction $D_k$ used by CauchyCos by multiplying on the right by the matrix ${(A)}^{-1}$, so it can be viewed as a right-preconditioning of the method CauchyCos. Notice also that ${\widehat{D}}_k$ resembles the residual direction $(X^{(k)}A-I)$ used in the minimal residual iterative method (MinRes) for minimizing ${\parallel I-XA\parallel }_F$ in the least-squares sense (see e.g., [6,10]). Nevertheless, the scaling factors in (7) differ from the scaling factors in the classical residual direction at $X^{(k)}$. Note that MinRes here should not be confused with the Krylov method MINRES [10].

For solving (3), we now present a variation of the CauchyCos Algorithm, that will be referred as the MinCos Algorithm, which from a given initial guess $X_0$ produces a sequence of iterates using the search direction ${\widehat{D}}_k$, while remaining in the compact set $S\cap T$. This new algorithm consists of simply replacing $-\nabla F(X^{(k)})$ in the CauchyCos Algorithm by ${\widehat{D}}_k$.

As before, we note that if we start from $X^{(0)}=(\sqrt{n}/{\parallel A\parallel }_F)I$ then by construction $\parallel X^{(k)}{A\parallel }_F=\sqrt{n}$, for all $k\ge 0$. For that initial guess, $trace(X^{(0)}A)=\langle X^{(0)}A,I\rangle >0$ and again by construction all the iterates remain in the $PSD$ cone. Notice also that, at each iteration, we now need to compute the two matrix–matrix products: $X^{(k)}A$, and ${\widehat{D}}_{k}A$, which for dense matrices require $n^3$ flops each. Every one of the remaining calculations (inner products and Frobenius norms) are obtained with n column-oriented inner products that require n flops each. Summing up, in the dense case, the computational cost of each iteration of the MinCos Algorithm is $2n^3+O(n^2)$ flops. In Section 2.5, we will discuss a sparse version of the MinCos Algorithm and its computational cost.

We start by noticing that, unless we are at the solution, the search direction ${\widehat{D}}_k$ is a descent direction.

We now establish the commutativity of all iterates with the matrix A.

It is worth noticing that using Lemma 2.9 and (5), it follows by simple calculations that $Z^{(k)}$, $X^{(k)}$, and $X^{(k)}A$ in the MinCos Algorithm are symmetric matrices for all k. These three sequences generated by the MinCos Algorithm are also uniformly bounded away from zero, and so the algorithm is well-defined.

Since $X^{(k)}\in S$, then $\sqrt{n}=\parallel X^{(k)}{A\parallel }_F\le \parallel X^{(k)}{A}^{1/2}{\parallel }_F{\parallel A^{1/2}\parallel }_F$, which combined with (8) implies that is bounded away from zero for all k. Moreover, since A is nonsingular then is also bounded away from zero for all k. ☐

For a given matrix A, the merit function $\Phi (X)=\frac{1}{2}{\parallel I-XA\parallel }_F^2$ has been widely used for computing approximate inverse preconditioners (see, e.g., [5,6,15,17,18,19,20,22]). In that case, the properties of the Frobenius norm permit in a natural way the use of parallel computing. Moreover, the minimization of $\Phi (X)$ can also be accomplished imposing a column-wise numerical dropping strategy leading to a sparse approximation of $A^{-1}$. Therefore, when possible, it is natural to compare the CauchyCos and the MinCos Algorithms applied to the angle-related merit function $F(X)$ with the optimal Cauchy method applied to $\Phi (X)$ (referred from now on as the CauchyFro method), and also to the Minimal Residual (MinRes) method applied to $\Phi (X)$ (see, e.g., [6,15]). Notice that, as we mentioned before in Section 2.3, with respect to CauchyCos and MinCos, MinRes can be seen as a right-preconditioning version of the method CauchyFro.

The gradient of $\Phi (X)$ is given by $\nabla \Phi (X)=-A^T(I-XA)$, and so the iterations of the CauchyFro method, from the same initial guess $X^{(0)}=(\sqrt{n}/{\parallel A\parallel }_F)I$ used by MinCos and CauchyCos, can be written as where $G_k=-\nabla \Phi (X^{(k)})$ and the step length ${\alpha }_k>0$ is obtained as the global minimizer of $\Phi (X^{(k)}+\alpha G_k)$ along the direction $G_k$, as follows where $R_k=I-A{X}^{(k)}$ is the residual matrix at $X^{(k)}$. The iterations of the MinRes method can be obtained replacing $G_k$ by the residual matrix $R_k$ in (9) and (10) (see [6] for details). We need to remark that in the dense case, the CauchyFro method needs to compute two matrix-matrix products per iteration, whereas the MinRes method by using the recursion $R_{k+1}=R_k-{\alpha }_{k}A{R}_k$ needs one matrix-matrix product per iteration.

We now discuss how to dynamically impose sparsity in the sequence of iterates $\left\{X^{(k)}\right\}$ generated by either the CauchyCos Algorithm or the MinCos Algorithm, to reduce their required storage and computational cost.

A possible way of accomplishing this task is to prescribe a sparsity pattern beforehand, which is usually related to the sparsity pattern of the original matrix A, and then impose it at every iteration (see e.g., [4,20,21,22]). At this point, we would like to mention that although there exist some special applications for which the involved matrices are large and dense [30,31], frequently in real applications the involved matrices are large and sparse. However, in general, the inverse of a sparse matrix is dense anyway. Moreover, with very few exceptions, it is not possible to know a priori the location of the large or the small entries of the inverse. Consequently, it is very difficult in general to prescribe a priori a nonzero sparsity pattern for the approximate inverse.

As a consequence, to force sparsity in our gradient related algorithms, we use instead a numerical dropping strategy to each column (or row) independently, using a threshold tolerance, combined with a fixed bound on the maximum number of nonzero elements to be kept at each column (or row) to limit the fill-in. This combined strategy will be fully described in our numerical results section.

In the CauchyCos and MinCos Algorithms, the dropping strategy must be applied to the matrix $Z_{k+1}$ right after it is obtained at Step 6, and before computing $X_{k+1}$ at Step 7. That way, $X_{k+1}$ will remain sparse at all iterations, and we guarantee that $X_{k+1}\in S\cap T$. The new Steps 7 and 8, in the sparse versions of both algorithms, are given by

7:

Apply numerical dropping to $Z^{(k+1)}$ with a maximum number of nonzero entries;

8:

Set $X^{(k+1)}=s\sqrt{n}\frac{Z^{(k+1)}}{\parallel Z^{(k+1)}{A\parallel }_F}$, where $s=1$ if $trace(Z^{(k+1)}A)>0$, $s=-1$ else.

Notice that, since all the involved matrices are symmetric, the matrix-matrix products required in both algorithms can be performed using sparse-sparse mode column-oriented inner products (see, e.g., [6]). The remaining calculations (inner products and Frobenius norms), required to obtain the step length, must be also computed using sparse-sparse mode. Using this approach, which takes advantage of the imposed sparsity, the computational cost and the required storage of both algorithms are drastically reduced. Moreover, using the column oriented approach both algorithms have a potential for parallelization.

To add understanding to the properties of the new CauchyCos and MinCos Algorithms, we start by testing their behavior, as well as the behavior of CauchyFro and MinRes, without imposing sparsity. Since the goal is to compute an approximation to $A^{-1}$, it is not necessary to carry on the iterations up to a very small tolerance parameter ϵ, and we choose $ϵ=0.01$ for our experiments. For all methods, we stop the iterations when $min\{F(X^{(k)}),\Phi (X^{(k)})\}\le ϵ$.

Table 2 shows the number of required iterations by the four considered algorithms when applied to some of the test functions, and for different values of n. No information in some of the entries of the table indicates that the corresponding method requires an excessive amount of iterations as compared with the MinRes and MinCos Algorithms. We can observe that CauchyFro and CauchyCos are not competitive with MinRes and MinCos, except for very few cases and for very small dimensions. Among the Cauchy-type methods, CauchyCos requires less iterations than CauchyFro, and in several cases the difference is significant. The MinCos and MinRes Algorithms were able to accomplish the required tolerance using a reasonable amount of iterations, except for the Lehmer(n) and minij(n) matrices for larger values of n, which are the most difficult ones in our list of test matrices. The MinCos Algorithm clearly outperforms the MinRes Algorithm, except for the Poisson 2D (n) and Poisson 3D (n) for which both methods require the same number of iterations. For the more difficult matrices and especially for larger values of n, MinCos reduces in the average the number of iterations with respect to MinRes by a factor of 4.

In Figure 1, we show the (semilog) convergence history for the four considered methods and for both merit functions: $F(X)$ and $\Phi (X)$, when applied to the Wathen matrix for $n=20$ and $ϵ=0.01$. Once again, we can observe that CauchyFro and CauchyCos are not competitive with MinRes and MinCos, and that MinCos outperforms MinRes. Moreover, we observe in this case that the function $F(X)$ is a better merit function than $\Phi (X)$ in the sense that it indicates with fewer iterations that a given iterate is sufficiently close to the inverse matrix. The same good behavior of the merit function $F(X)$ has been observed in all our experiments.

Based on these preliminary results, we will only report the behavior of MinRes and MinCos for the forthcoming numerical experiments.

We now build sparse approximations by applying the dropping strategy, described in Section 2.5, which is based on a threshold tolerance with a limited fill-in ($lfil$) on the matrix $Z^{(k+1)}$, at each iteration, right before the scaling step to guarantee that the iterate $X^{(k+1)}\in S\cap T$. We define $thr$ as the percentage of coefficients less than the maximum value of the modulus of all the coefficients in a column. To be precise, for each i-th column, we select at most $lfil$ off-diagonal coefficients among the ones that are larger in magnitude than $thr\times \parallel {(Z^{(k+1)})}_i{\parallel }_{\infty }$, where ${(Z^{(k+1)})}_i$ represents the i-th column of $Z^{(k+1)}$. Once the sparsity has been imposed at each column and a sparse matrix is obtained, say $sp(Z^{(k+1)})$, we guarantee symmetry by setting $Z^{(k+1)}=(sp(Z^{(k+1)})+sp{(Z^{(k+1)})}^T)/2$.