# Matlab help related with the scaled Newton's iteration method

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.

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