# Convergence of iterative algorithm.

For quite a long time I'm trying to prove convergence of an iterative algorithm in case of a particular system of nonlinear equations.

Here are some characteristics of this system:
It consists of n equation and n variables.
Every equation is in similar form - sum of products = constant.
The lenght of every product is the same (it will be denoted as k).
The number of elements in sum might be different in each equation.
In every product in equation $i$ one of elements is $x_{i}$ - this is very important.
This system is "symmetrical", it means that if $x_{i} \cdot x_{j} \cdot ...$ is one of elements of equation i, then it is also in equaton j.
$b_{i} > 0$ - where $b_{i}$ is intercept in equation i.

I'll write an example of such equation for k=3 and n=6:
$x_{1} \cdot x_{3} \cdot x_{6} + x_{1} \cdot x_{2} \cdot x_{4} = b_{1}$
$x_{2} \cdot x_{1} \cdot x_{4} + x_{2} \cdot x_{5} \cdot x_{6} = b_{2}$
$x_{3} \cdot x_{1} \cdot x_{6} + x_{3} \cdot x_{4} \cdot x_{6} = b_{3}$
$x_{4} \cdot x_{1} \cdot x_{2} + x_{4} \cdot x_{3} \cdot x_{6} = b_{4}$
$x_{5} \cdot x_{2} \cdot x_{6} = b_{5}$
$x_{6} \cdot x_{1} \cdot x_{3} + x_{6} \cdot x_{2} \cdot x_{5} + x_{6} \cdot x_{3} \cdot x_{4} = b_{6}$

It is very easy to transform this equation to following form (it just need to be divided once).
$x_{i} = b_{i} / something$ , $x_{i}$ is only on left-hand side of ith equation.

If we have all equation in such form then the fixed point is solution of it.
I've experimentally checked that algorithm analogic to Gauss-Seidel is covergent (i've checked ~100 random examples, and in every case it was convergent).
By analogic to Gauss-Seidel algorithm I mean:
1) Choose any initial solution $[x_{1}^{0} , ... , x_{n}^{0}]$
2.1) Calculate value of $x_{1}^{i+1}$ using $[x_{2}^{i} , ... , x_{n}^{i}]$
2.2) Calculate value of $x_{2}^{i+1}$ using $[x_{1}^{i+1} , ... , x_{n}^{i}]$
...
2.n) Calculate value of $x_{n}^{i+1}$ using $[x_{1}^{i+1} , ... , x_{n-1}^{i+1}]$
3) If solution is good enough stop, otherwise go to 2.1

I've tried Banach fixed point theorem, but is hard to say anything about spectral radius. Does anyone have a clue how to prove convergence of this algorithm?

Edited 7.02.2009 11:09
I've found another restriction. If we denote by $m$ number of all products (in this example it would be 12, becuase in 1st,2nd,3th,4th we have 2 products, in 5th we have 1, and in 6th we have 3). Then following equation is true:
$m = k \cdot \sum_{i=1}^{n}{b_{i}}$ Which also implies thah $\sum_{i=1}^{n}{b_{i}}$ is a natural number (from the symmetry), but it doesn't imply that any of $b_{i}$ is natural.

-
Why tag this as nt? – lhf Feb 4 '10 at 17:53
Sorry it was mistake, I've clicked on number theory instead of numerical analysis. – Tomek Tarczynski Feb 4 '10 at 18:02
I bet you wanted all $x_j>0$ everywhere, but since you have never said it explicitly, could you please confirm? Also the system $x_1x_2=b_1$, $x_2x_1=b_2$ satisfies all your conditions but I wish you good luck solving it iteratively or otherwise with $b_1\ne b_2$. Are we just assuming that the solution exists? – fedja Feb 4 '10 at 22:41
Both of your assumptions arue true. - $x_{j} > 0$ - And we assume that solution exists. – Tomek Tarczynski Feb 5 '10 at 8:39

One idea we'll certainly need here is that the sum of absolute values of the discrepancies $D_k$ (left hand side minus right hand side) does not increase after each operation due to the symmetry property. If you move $x_1$ by $\delta$ towards the solution of the first equation, whatever shift you have in the $k$th equation because of that, it also arises in the first equation with the sign that decreases the absolute value of the first discrepancy. Now, the problem is that it may never help because in every other equation the sign of the discrepancy is opposite to the sign in the first equation, so whatever you gain in equation 1, you lose the same amount in other equations combined. This is exactly what happens in my trivial counterexample. Note that in this case, if we originally have $D_1>0$, then all other $D_k\le 0$ (otherwise we gained somewhere instead of loosing). But then when we move $x_1$ down to kill $D_1$, we move every other $D_k$ to the negative domain whence at the second iteration we'll gain in $D_3$ (if the equations are more than 2 and the 2-system is unique).

The real question is why you cannot escape to infinity instead of converging to a solution. Note that the convergence of discrepancies to $0$ is always geometric (we gain at least some multiple of the maximal initial discrepancy when we go over the full cycle), so your algorithm always stops reasonably fast and outputs something even if there is no solution making your random system experimental check totally useless. Still, if the iterations do not converge, we must have the smallest variable to go to to 0 and the largest one to infinity geometrically as well. That is certainly possible: you clearly have each LHS dominated by the sum of the others and can choose the RHS so that one number is greater than the sum of all others, so there are equations of every order without solutions and if we had stayed in the compact domain (or even returned to it infinitely many times), every limit point would be a solution due to the above argument. Together with the fact that you cycle in the 2 system even when the right hand sides are equal, this makes me not so sure that mere existence of a solution is a strong enough assumption. It even seems quite feasible that you may have a solution with a small attracting basin and everything else just goes away.

Now, are you doing it just for fun, or you really need it for something?

-
I'm still analysing your answer, but I can definietly say that it's for something - all variables has economical interpretation. – Tomek Tarczynski Feb 6 '10 at 9:23
It is actually nonsense when $k > 2$, but when $k=2$ and each variable occurs in at least 2 products, what I said is correct (with minor modifications), so I'll keep it here. Let's try to understand the k=2 case fully first. – fedja Feb 6 '10 at 14:28
Yesterday I've been trying to prove modified version of this algorithm. Instead of doing this steps consecutive for $x_{1}, x_{2} ,..$ I was always doing them for $x_{i}$ with highest asbolute value of $D_{i}$ (the convergance is even faster in this case). From the empirical point of view I may say that after doing step for $x_{i}$ sum of 'new' $D_{i}$ might be higher than the provious one, I was checking it when k=2. – Tomek Tarczynski Feb 7 '10 at 10:05

Your system of equations is of the form $Az = b$, where $z_{i(\alpha)} = x^\alpha$ and $i(\alpha)$ is the grlex index of

$\alpha \in$ $X_{n,k} \equiv$ {$\beta \in \mathbb{Z}^n: \sum_j \beta_j = k$}.

(See Uniquely generate all permutations of three digits that sum to a particular value? for information on this index.)

Solving $Az = b$ is tantamount to solving the equations (if a solution exists), because it is trivial to obtain $x$ from $z$. So really you just need to solve a matrix equation once you sort through the indexing issues.

-
But the linear system above is underdetermined, is it not? And then, the problem of finding x from z is overdetermined, or so it seems to me. So this is not a great way to decompose the problem. – Harald Hanche-Olsen Feb 4 '10 at 19:37
Fair enough, I wasn't thinking all the way through. But using this index I can do the following in MATLAB... A = zeros(56); A(1,42) = 1; A(1,48) = 1; A(2,48) = 1; A(2,22) = 1; A(3,42) = 1; A(3,14) = 1; A(4,48) = 1; A(4,14) = 1; A(5,22) = 1; A(6,42) = 1; A(6,22) = 1; A(6,14) = 1; ...and then pA=pinv(A);sparse(pA.*(abs(pA)>10^-10)) gives a pretty tractable psuedoinverse. – Steve Huntsman Feb 4 '10 at 20:30
BTW, the commands above are related to the specific example in the question. – Steve Huntsman Feb 4 '10 at 20:30
However, the rank of A in this example is 4 rather than 6, so there are certainly some uniqueness issues. – Steve Huntsman Feb 4 '10 at 20:35
Please tell me matlab has a less verbose way of entering sparse matrices! :) – Mariano Suárez-Alvarez Feb 4 '10 at 20:36