Matlab help related with the scaled Newton's iteration method

Question

I need a little help with Matlab code for the method mentioned in this paper for computing the inverse of the matrix.

In this paper, the author presented scaled Newton iteration given by $X_{k+1} = \alpha_{k+1}(2I - X_{k}A)X_{k}, k = 0, 1, 2, \ldots$ for a given matrix $A$.

Initial approximation $X_{0} = \alpha_{0}A^t$ , $\alpha_{0} = \frac{2}{\sigma_{\min}+\sigma_{\max}}$ , $\sigma_{\min}$, $\sigma_{\max}$ are the the minimum and maximum singular values of the matrix $A$.

$\alpha_{k+1} = \frac{2}{1+(2 - \rho_{k} )\rho_{k}}$

$\rho_{0} = \alpha_{0} \sigma_{\min}$

$\rho_{k+1} = \alpha_{k+1} (2 - \rho _{k} )\rho_{k} $

Here is the algorithm that I want to implement in Matlab.

1. Find singular values $\sigma_1 \geq \sigma_2\ldots \geq\sigma_{r}$ of $A$
2. $\alpha_{0} = \frac{2}{\sigma_{r}^{2}+\sigma_{1}^2}$
3. $\rho_{0} = \alpha_{0} \sigma_{r}$
4. $X_{0} = \alpha_{0}A^t$
5. for $k = 1, 2\ldots\quad$: $\alpha_k = \frac{2}{[1 + (2 - \rho_{k-1})\rho_{k-1}]}$
6. $X_{k} = \alpha_{k}(2I - X_{k-1}A)X_{k-1}$
7. $\rho_{k} = \alpha_{k} (2 - \rho_{k-1})\rho_{k-1}$

Here is the matlab code that I have made

A = rand(5);   % given matrix 
  svds = svd(A)  % singular value of the matrix A
   s = min(svds)  %smallest singular value
   s1  = max(svds) %lsrgest singular value
   s2 = s^2 + s1^2  
   alpha0 = 2/s2  %alphazero
  rho0 = alpha0 *s
  x0 = alpha0*A';  %initial approximation

  I = eye(5)

  iter = 0
   for i = 1:20
       a1 = 2/(1 + (2 - rho0)*rho0)
       g1 = x0+x0*(I - A*x0);
       x1 = a1*(g1);     % approximation of inverse
       rho1 = a1*(2*rh0 -rh0*rh0 )

       x0 = g1;
       rho0 = rho1

       iter = iter +1
   end 

This code gives the correct result. But after few iterations, the values of rho1 and a1 coincide and become 1 and, hence, it is taking an equal number of iterations compared to the classical Newton iteration. While it should take less number of iterations as mentioned in the paper. I want to know where I am wrong.

Answer 1

I can't access the article, but from what you describe, it looks quite normal that $\alpha$ and $\rho$ become close to $1$: if the method converges, you will have

$$ X = \alpha(2I - XA)X $$

which is equivalent to (supposing $X$ invertible)

$$ I = \alpha(2I - XA) $$

so

$$ XA = (2-\frac{1}{\alpha}) I$$

and (supposing $A$ is invertible)

$$ X = (2-\frac{1}{\alpha})A^{-1}$$

If you want $X$ to converge to $A^{-1}$, you need $\alpha = 1$ and since $\alpha = \frac{2}{1+(2-\rho)\rho}$, you need $\rho$ to go to $1$ as well.

To me, it looks like the method is just working as it should. The reduced number of iterations might come from the "warm-up" phase at the beginning (that's just a guess).