The most obvious way to solve a linear equation is to apply the Gauss
elimination method [PFTV88]. Unfortunately
it fails to solve the radiosity equation for more complex models
effectively, since it has 
complexity, and also it
accumulates the round of errors of digital computers and magnifies these
errors to the extent that the matrix is close to singular.
Fortunately another technique, called iteration, can overcome both problems. Examining the radiosity equation,
we will see that it gives the equality of the energy
which has to be radiated due to emission and reflection (right side) and the
energy really emitted (left side). Suppose that only estimates are available
for 
radiosities, not exact values. These estimates can be
regarded as right side values, thus having substituted them into the
radiosity equation, better estimates can be expected on the left sides. If
these estimates were exact -- that is they satisfied the radiosity
equation --, then the iteration would not alter the radiosity values. Thus,
if this iteration converges, its limit will be the solution of the original
radiosity equation.
In order to examine the method formally, the matrix version of the radiosity equation is used to describe a single step of the iteration:
A similar equation holds for the previous iteration too. Subtracting the two equations, and applying the same consideration recursively, we get:
The iteration converges if
for some matrix norm. Let us use the 
norm defined as the maximum of absolute row sums
and a vector norm that is compatible with it:
Denoting 
by q, we have:
according to the properties of matrix norms. Since 
represents the
portion of the radiated energy of surface i, which actually reaches
surface j, 
is that portion which is
radiated towards any other surface. This obviously cannot exceed 1, and for
physically correct models, diffuse reflectance 
, giving a
norm that is definitely less than 1. Consequently q<1, which
provides the convergence with, at least, the speed of a geometric series.
The complexity of the iteration solution depends on the operations
needed for a single step and the number of iterations providing convergence.
A single step of the iteration requires the multiplication of an N
dimensional vector and an 
dimensional matrix, which
requires 
operations.
Concerning the number of necessary steps, we concluded that the speed of
the convergence is at least geometric by a factor 
. The
infinite norm of 
is close to being
independent of the number of surface elements, since as the number of
surface elements increases, the value of form factors decreases,
sustaining a constant sum of rows, representing that portion of the energy
radiated by surface i, which is gathered by other surfaces,
multiplied by the diffuse coefficient of surface i. Consequently, the
number of necessary iterations is independent of the number of surface
elements, making the iteration solution an process.