Quickly and efficiently calculating an eigenvector

Short version of my question:

What would be the optimal way of calculating an eigenvector for a matrix A, if we already know the eigenvalue belonging to the eigenvector?

Longer explanation:

I have a large stochastic matrix A which, because it is stochastic, has a non-negative left eigenvector x (such that A^Tx=x).

I'm looking for quick and efficient methods of numerically calculating this vector. (Preferrably in MATLAB or numpy/scipy - since both of these wrap around ARPACK/LAPACK, any one would be fine).

I know that 1 is the largest eigenvalue of A, so I know that calling something like this Python code:

from scipy.sparse.linalg import eigs
vals, vecs = eigs(A, k=1)

will result in vals = 1 and vecs equalling the vector I need.

However, the thing that bothers me here is that calculating eigenvalues is, in general, a more difficult operation than solving a linear system, and, in general, if a matrix M has eigenvalue l, then finding the appropriate eigenvector is a matter of solving the equation (M - 1 * I) * x = 0, which is, in theory at least, an operation that is simpler than calculating an eigenvalue, since we are only solving a linear system, more specifically, finding the nullspace of a matrix.

However, I find that all methods of nullspace calculation in MATLAB rely on svd calculation, a process I cannot afford to perform on a matrix of my size. I also cannot call solvers on the linear equation, because they all only find one solution, and that solution is 0 (which, yes, is a solution, but not the one I need).

Is there any way to avoid calls to eigs-like function to solve my problem more quickly than by calculating the largest eigenvalue and accompanying eigenvector?

Here's one approach:

Let x denote the (row) left^† eigenvector associated to eigenvalue 1. It satisfies the system of linear equations (or matrix equation) xA = x, or x(A−I)=0.
To avoid the all-zeros solution to that system of equations, remove the first equation and arbitrarily set the first entry of x to 1 in the remaining equations.
Solve those remaining equations (with x₁ = 1) to obtain the other entries of x.

Example using Matlab:

>> A = [.6 .1 .3
        .2 .7 .1
        .5 .1 .4]; %// example stochastic matrix
>> x = [1, -A(1, 2:end)/(A(2:end, 2:end)-eye(size(A,1)-1))]
x =
   1.000000000000000   0.529411764705882   0.588235294117647
>> x*A %// check
ans =
   1.000000000000000   0.529411764705882   0.588235294117647

Note that the code -A(1, 2:end)/(A(2:end, 2:end)-eye(size(A,1)-1)) is step 3.

In your formulation you define x to be a (column) right eigenvector of A^T (such that A^Tx = x). This is just x.' from the above code:

>> x = x.'
x =
   1.000000000000000
   0.529411764705882
   0.588235294117647
>> A.'*x %// check
ans =
   1.000000000000000
   0.529411764705882
   0.588235294117647

You can of course normalize the eigenvector to sum 1:

>> x = x/sum(x)
x =
   0.472222222222222
   0.250000000000000
   0.277777777777778
>> A.'*x %'// check
ans =
   0.472222222222222
   0.250000000000000
   0.277777777777778

^{^† Following the usual convention. Equivalently, this corresponds to a right eigenvector of the transposed matrix; see below.}