Computing/certifying bounds for eigenvalue computations for big sparse 0-1 matrices I am finishing my PhD, and at some point I give some bounds on let's say a constant $\lambda$. This constant is the maximum eigenvalue of a matrix $M$ with only zeroes and ones (so the Perron-Frobenius theorem applies and says that there is a simple maximum eigenvalue $\lambda$).
To compute it, I choose some vector $V$ (which has a meaning in my problem, it cannot be orthogonal to the eigenvector for $\lambda$). Then I computes the ratios $ \|M^nV\|/\|M^{n-1}V\|$. This should converge towards $\lambda$.
I do around $n=100$ iterations, and it converges already from ~30. However, I have no proof that it has begun converging, that my approximation is close to $\lambda$. I would like to certify that some decimals are the true ones, i.e. that my approximation is good up to $1/k$ for some value of $k$ I choose.
My matrices can have ~1 million lines, but are sparse, hence I am able to compute matrix-vector products for the "power method" I am using to compute the maximum eigenvalue.
Would anyone know how to do that, and refer to some paper proving it?
Another approach would be, since I obtain a vector $V_0$ should not be far from an eigenvector for $\lambda$, to certify that it is not far from the actual eigenvector (because we can see that $MV_0 \approx \lambda V_0$) and that the value we compute is "hence" (to be proved) not far from $\lambda$.
Thanks a lot in advance,
Alexandre
PS: it has no importance here, but just for the background, I am computing bounds on the growth rates of some (various types of) dominating sets in rectangular grids.
 A: PART 1. 
Firstly, your matrix is only non-negative; consequently, $\lambda=\rho(A)$ is an eigenvalue of $A$ associated to a non-negative eigenvector $u$; yet, $\rho(A)$ is not necessarily simple and there may be other eigenvalues ​​of A of same modulus. See for example $A=\begin{pmatrix}0&0&1\\1&0&0\\0&1&0\end{pmatrix}$.  Anyway, the first thing to do is to randomly choose the vector $V$ among the NON-NEGATIVE vectors.
Here your matrices are very large; then, with a great probability, $\lambda$ is simple $>0$ and is the unique eigenvalue with modulus $\rho(A)$. Now, we consider the second largest eigenvalue (in modulus); a random real polynomial has much more complex roots than real ones (the quotient is in $O(\log(n)/n)$; then, with a great probability, the second largest eigenvalues are not real and are, generically, a pair ($Av=\mu v,A\overline{v}=\overline{\mu}\overline{v})$ (we assume that we are in this case but, possibly, $\mu$ may be real).
Let $V=au+bv+\overline{b}\overline{v}+\cdots$ be the decomposition of $V$ on a basis $u,v,\overline{v},\cdots$; notice that $a,b$ are non-zero with probability $1$. Since $M^kV\approx \lambda^k(au+2Re((\mu/\lambda)^kbv))$ and $\lambda>0$, we can deduce that $M^kV/||M^kV||$ tends (in the sphere) to an eigenvector of $A$. On the other hand,
$\tau_k=||M^{k+1}V||/||M^kV||\approx \lambda||au+2Re((\mu/\lambda)^{k+1}bv)||/||au+2Re((\mu/\lambda)^{k}bv)||$
$\approx \lambda(1+O_1((\mu/\lambda)^{k})/(1+O_2((\mu/\lambda)^{k})\approx \lambda(1+O(|\mu/\lambda|^{k})$.
Finally, $\tau_k$ tends to $\lambda$ and the approximation of $\lambda$ by $\tau_k$ has (roughly speaking) $l$ significative digits where $l$ is given by the relation $|\mu/\lambda|^k=10^{-l}$.
Notice that $v,\overline{v}$ generate a real plane $\Pi$ that is $A$-invariant. To find an estimate of the error about $(\tau_k,V)$, it suffices to obtain approximations -with few digits- of a basis of $\Pi$.
REMARKS.  i) An iteration must be written as follows
$W:=AV:V:=\dfrac{1}{||W||}W$:
ii) A good new. Generically, for large random matrices,  $|\mu/\lambda|$ is small ($< 1/10$ for $n\geq100$). This is why convergence is done in $\approx 30$ iterations.
iii) Yet, you want a certified stop. The bad new is that the third eigenvalue-block (in modulus) may be close to the second one. You can proceed as follows
PART 2.
Calculate an approximation of "the" eigenvector $V_1$ of $A^T$ associated to $\lambda$. As above, an iteration is
$W:=A^TV_1:V_1:=\dfrac{1}{||W||}W$.
In a first step, we can stop the iterations when the first $5$ digits do not change during $5$ successive iterations (for example).
Then $V_1^{\perp}$ is $A$-invariant. Randomly choose $V_2\in V_1^{\perp}$ and write the iteration (we force it to stay in $V_1^{\perp}$)
$W:=AV_2:V_2:=\dfrac{1}{||W||}W:V_2:=V_2-\dfrac{V_2^TV_1}{||V_1||^2}V_1$.
Thus $\{V_2,V_3=AV_2\}$ is an approximation of a basis of $\Pi$. According to Remark 3, the convergence may be slow; to obtain $3$ significand digits, it may be necessary to perform at least $200$ iterations...
Calculate, by least squares, $a,b$ s.t. $AV_3\approx aV_2+bV_3$. Then the matrix of $A_{|\Pi}$ is $Z=\begin{pmatrix}0&a\\1&b\end{pmatrix}$.
Finally, the relative error is $\dfrac{1}{\lambda^k}||Z^k||$. It remains to solve $\dfrac{1}{\lambda^k}||Z^k||=10^{-l}$.
