In this paper, we present a novel algorithmic approach, the hybrid matrix geometric/invariant subspace method, for finding the stationary probability distribution of the finite quasi-birth-death (QBD) process which arises in performance analysis of computer and communication systems. Assuming that the QBD state space is defined in two dimensions with m phases and K+1 levels, the solution vector for level k, /spl pi//sub k/, 0/spl les/k/spl les/K is shown to be in a modified matrix geometric form /spl pi//sub k/=/spl upsi//sub 1/R/sub 1//sup k/+/spl upsi//sub 2/R/sub 2//sup K-k/ where R/sub 1/ and R/sub 2/ are certain solutions to two nonlinear matrix equations and /spl upsi//sub 1/ and /spl upsi//sub 2/ are vectors to be determined using the boundary conditions. We show that the matrix geometric factors R/sub 1/ and R/sub 2/ can simultaneously be obtained independently of K via finding the sign function of a real matrix by an iterative algorithm with quadratic convergence rates. The time complexity of obtaining the coefficient vectors /spl upsi//sub 1/ and /spl upsi//sub 2/ is shown to be O(m/sup 3/ log/sub 2/ K) which indicates that the contribution of the number of levels on the overall algorithm is minimal. Besides the numerical efficiency, the proposed method is numerically stable and in the limiting case of K/spl rarr//spl infin/, it is shown to yield the well-known matrix geometric solution /spl pi//sub k/=/spl pi//sub 0/R/sub 1//sup k/ for infinite QBD chain.