Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I simply would like to find the inverse of a given matrix $A$ by the iterative method $X_{k+1}=X_k(2I-AX_k)$ using $X_0=\frac{1}{\sigma_{max}^2}A^*$. An example is as follows

SeedRandom[123]; m = n = 500;
A = RandomReal[10, {m, n}]; Id = SparseArray[{{i_, i_} -> 1.}, {m, m}];
tolerance = 10^-6; max = 50; k = 0; R[0] = 1;
X[0] = (1./SingularValueList[A, 1][[1]]^2)*Transpose[A];
While[k < max && R[k] >= tolerance,
  XX = A.X[k]; X[k + 1] = X[k].(2 Id - XX);
  R[k + 1] = Norm[Id - X[k + 1].A, Infinity];
Table[R[l], {l, 1, k}]

As can be obsereved the number iterations is too high for a $500\times500$ matrix. What I want to do now is to scale this method by introducing the parameter $\alpha_k$ per cycle, i.e. $X_k<=\alpha_kX_k$. In fact, $\alpha_k$ must be computed adaptively so as to solve the following optimization problem $\|I-\alpha_kAX_k\|_2<\|I-AX_k\|_2$. In this way, we might be able to reduce the number of iterations. Note that $\alpha_k$ must rapidly tend to 1 in order not to cause instability.

I will be thankful if anyone lend me a hand in caluculating the scaling parameter per cycle.

share|improve this question
+1 for an interesting question, although I'm not sure it's on topic exactly. You might find you receive more helpful suggestions at Computational Science.SE, which specializes in this type of question. – Oleksandr R. Jun 13 '14 at 0:34
For searching purposes: the iterative method being discussed here is referred to as Newton-Schulz iteration. – J. M. May 27 '15 at 11:19

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.