3
$\begingroup$

I want $S^k$, with $S=I-\Lambda^{-1}M$, to tend to zero quite fast as $k\rightarrow \infty$, as this is what drives the convergence in a fixed-point algorithm. Here $M=X^TX$ is a fixed $m\times m$ matrix, $I$ is the $m\times m$ identity matrix, $\Lambda$ is an $m\times m$ diagonal matrix, and $X$ is an $n \times m$ matrix so that $M=X^TX$ is $m\times m$, symmetric and positive semidefinite. I am trying to find a good $\Lambda$ that achieves this goal, yet one that is simple easy to compute.

I know the convergence speed is driven by the largest eigenvalue of $S$. That eigenvalue must be $<1$ in absolute value. Let $\lambda_i$ be the $i$-th element of $\Lambda$. If $\Lambda$ is chosen so that the resulting elements of $S$ are in some sense, close enough to zero - as "close" as they can be - then one would expect fast convergence, and it does work in practice. For $\lambda_i$, I chose the diagonal element of $M$ in the $i$-row, divided by the sum of the squares of the elements of $M$ in the $i$-th row. I am moderately happy with the results (at least for matrices up to $m=6$) but I am wondering if it is possible to get better $\lambda_i$'s, that on average will further boost convergence to zero. Also, I want to keep the $\lambda_i$'s as simple as possible. All elements in all the matrices are real numbers.

My question: Does my choice of $\Lambda$ always lead to $S^k\rightarrow 0$ ($k\rightarrow\infty$), and is there a better choice (yet as simple as possible) that will make convergence faster?

$\endgroup$

2 Answers 2

2
$\begingroup$
  1. As Iosif Pinelis pointed out, you presumably meant $$ \lambda_i^{-1} = \frac{e_i^T M e_i}{\lVert M e_i \rVert_2^2}, $$ where $e_i$ is the $i$th standard basis vector, which is indeed a reasonable choice, and can be motivated as being the unique minimizer of the Frobenius norm $$ \lVert I - \Lambda^{-1} M \rVert_F $$ under the constraint that $\Lambda$ be diagonal. You can find some literature on preconditioners found by minimizing this sort of Frobenius norm subject to some sparsity constraint under the keywords "sparse approximate inverse" (SAI). The specific choice above, where the preconditioner is constrained to be diagonal, was considered in [1] under the name "SPAI-0". In the context of [1], this preconditioner was recommended as a component (the smoother) of a more complicated preconditioner (multigrid), precisely because for their application the spectral radius of $S$ was expected to still be close to 1, implying slow convergence.

  2. In more generality, it is a standard result in iterative methods for solving linear systems that the optimal choice of damping parameter $\omega$ minimizing the spectral radius of $$ I - \omega B A $$ for symmetric positive definite $A$ and $B$ is given by $$ \omega = \frac{2}{\lambda_{\text{min}} + \lambda_{\text{max}}}, $$ where $\lambda_{\text{min}}$ and $\lambda_{\text{max}}$ are the smallest and largest eigenvalues of $BA$, which are necessarily all positive real. This choice gives the spectral radius $$ \rho(I - \omega B A) = \frac{\kappa - 1}{\kappa + 1} < 1, \quad\text{where } \kappa = \frac{\lambda_{\text{max}}}{\lambda_{\text{min}}} . $$ Relating this back to your problem, any choice of SPD $\Lambda$, not just diagonal, can be made to converge with a suitable scaling, provided $M$ is nonsingular. This generalizes part 4 of Iosif's answer, which corresponds to the choice $B=I$ above. Note this argument reduces the problem to finding an optimal preconditioner $B$ ($\Lambda^{-1}$ in your case), in the sense of minimizing the condition number $\kappa$.

  3. It is easy to analyze the $2 \times 2$ case in full detail. Consider the matrix $$ M = \begin{bmatrix} 1 & \tfrac{\kappa-1}{\kappa+1} \\ \tfrac{\kappa-1}{\kappa+1} & 1 \end{bmatrix}, $$ which has eigenvalues $2/(\kappa+1)$ and $2\kappa/(\kappa+1)$, and condition number $\kappa$. One can show that the spectral radius of $S$ is minimized in this case by $\Lambda = I$ (assuming diagonal $\Lambda$, of course), and the minimum spectral radius is, again, $$ \rho(S) = \frac{\kappa - 1}{\kappa + 1} . $$ This shows, as Iosif pointed out, that the convergence can be arbitrarily slow for poorly conditioned $M$, even for the optimal choice of diagonal preconditioner, already in the $2 \times 2$ case. This also shows that the nonconvergence that occurs when $M$ is singular can be approached "continuously" in the limit $\kappa \to \infty$, in the sense that $\rho(S) \to 1$ continuously as $\kappa \to \infty$. Note that your choice of $\Lambda$, the one minimizing $\lVert I - \Lambda^{-1} M \rVert_F$, is suboptimal here, though it still gives $\rho(S) < 1$. In the worst case, it gives a $\rho(S)$ about 21% higher than optimal for this particular $M$.

[1] Bröker, Oliver; Grote, Marcus J., Sparse approximate inverse smoothers for geometric and algebraic multigrid, Appl. Numer. Math. 41, No. 1, 61-80 (2002). ZBL0995.65129.

$\endgroup$
2
  • $\begingroup$ This is useful information, with good generalizations. $\endgroup$ Apr 29, 2022 at 20:10
  • $\begingroup$ Great answer! Many thanks. $\endgroup$ Apr 29, 2022 at 21:33
2
$\begingroup$

$\newcommand{\La}{\Lambda}\newcommand{\la}{\lambda}$

  1. No, your choice of $\La$ will not always lead to $S^k\to0$. In fact, for some symmetric positive-semidefinite matrices $M$, no choice of $\La$ will lead to $S^k\to0$. In particular, suppose that $M$ is singular -- which will always be the case when $m>n$. Then the matrix $\La^{-1}M$ is also singular and hence has eigenvalue $0$, so that $S=I-\La^{-1}M$ will have eigenvalue $1$, which will preclude $S^k\to0$.

  2. Your choice of $\La$ seems strange. E.g., suppose that $M=2I$. Then, of course, one can choose $\La=2I$, to make $S=0$. On the other hand, here your $S$ will be $-3I$, so that your $S^k$ will exponentially diverge.

  3. To make $S=I-\La^{-1}M$ small, it seems reasonable to make the diagonal matrix $\La$ as close to $M$ as possible, in some sense. A simplest choice of this kind is to let the diagonal entries of $\La$ be the same as those of $M$. Then, if e.g. the diagonal entries of $M$ are greater than the sums of the non-diagonal entries of the corresponding rows, then the $\ell^\infty\to\ell^\infty$ operator norm of $S$ will be $<1$, so that we will indeed have $S^k\to0$.

  4. On a much more positive note, if the symmetric positive-semidefinite matrix $M$ is non-singular (and hence positive definite), then there always is a choice of $\La$ such that $S^k\to0$.

Indeed, let $0<\la_1\le\cdots\le\la_m$ be the eigenvalues of the symmetric positive-definite matrix $M$. Let \begin{equation} T:=I-U,\quad U:=\La^{-1/2}M\La^{-1/2},\quad\La:=\frac{\la_1+\la_n}2\,I. \end{equation} Then for the spectral/operator norm $\|U\|$ of $U$ we have \begin{equation} \|U\|\le\|\La^{-1/2}\|\,\|M\|\,\|\La^{-1/2}\|=\frac{2\la_n}{\la_1+\la_n}<2. \end{equation} Since the matrix $U$ is symmetric positive definite, it follows that $0<U<2I$, where $A<B$ means that $B-A$ is a symmetric positive definite matrix. So, \begin{equation} -I<T<I, \end{equation} which implies $T^k\to0$. Finally, since $S=\La^{-1/2}T\La^{1/2}$, we have \begin{equation} S^k=\La^{-1/2}T^k\La^{1/2}\to0, \end{equation} as claimed.

$\endgroup$
6
  • $\begingroup$ Thank you. I will add context to my question. In short I am solving a linear regression problem with fixed point algorithm; the goal is to solve other related problems in the future with this approach. Here $X$ is the data, so usually $n>m$. And if $M$ is singular, linear regression has issues too. This happens when one dependent variable is a linear combination of the other dependent variables. $\endgroup$ Apr 29, 2022 at 4:33
  • $\begingroup$ (continued) Barring these situations, not sure if my choice of $\Lambda$ always leads to convergence. It did work in all the examples tested, though sometimes not as nicely as expected. I tried your #3, but it does not work as well in general. My choice of $\Lambda$ is to minimize the sum of the squares of the coefficients in each row of $S$. $\endgroup$ Apr 29, 2022 at 4:36
  • $\begingroup$ It seems that for $1/\lambda_i$ (not $\lambda_i$), you chose the diagonal element of $M$ in the $i$-row, divided by the sum of the squares of the elements of $M$ in the $i$-th row. $\endgroup$ Apr 29, 2022 at 4:45
  • $\begingroup$ Yes, this is correct. Not saying it is a great choice. Sometimes the largest eigenvalue in absolute value is close to 1 (sometimes very close but always below 1 unless my computation is wrong) but even then, after 15 iterations of the fixed point, I get predictors with R-squared just as good as the exact solution blended with little noise (on synthetic data). I tried an adaptive $\Lambda$ (updated at each iteration) but so far it works 80% of the time and fails miserably on many examples. $\endgroup$ Apr 29, 2022 at 5:13
  • $\begingroup$ Take a look at just added item 4. I think this should work better. On the other hand, $M$ can be ill-conditioned (as likely to happen for large $m$), and then hardly anything will work well, I think. $\endgroup$ Apr 29, 2022 at 5:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.