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

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.


1 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) $$


$$ 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).


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.