This is the basis of a pretty common set of techniques to find ground state properties. The hard part is writing down the matrix and multiplying it against trial wavefunctions in a large many-body basis. The projection intuition itself is not enough, but it turns out we can use:
Projector Quantum Monte Carlo (there's a lot of literature on this, but see for example http://arxiv.org/abs/0807.0682 section IV) to sample efficiently the effect of hitting a trial state with high powers of the hamiltonian matrix.
Time-evolved Block Decimation in imaginary time. This technique is closely related to DMRG, and again it's just a matter of having a good trial state (a Matrix Product State) and an efficient way of applying the thermal evolution operator (technical, but all the details are in http://arxiv.org/abs/quant-ph/0301063)
I believe that it is true as long as there does not exist a non-trivial unitary operator $U$ that commutes with the Hamiltonian ($[H,U] = 0$) in the subspace of ground states. If such an operator exists then for a ground state $|\phi_0\rangle$ with energy $E_0$ we have $$HU|\phi_0\rangle = UH|\phi_0\rangle = E_0\left(U|\phi_0\rangle\right)$$ and so $U|\phi_0\rangle$ also has the lowest possible energy $E_0$ and it thus also a ground state. Note that the statement of non-triviality of $U$ is important. It needs to be non-trivial in the subspace of ground states, that is $U|\phi_0\rangle \neq e^{i\theta}|\phi_0\rangle$ for any phase $\theta$, otherwise there is no degeneracy. (Unitarity is needed so that $U|\phi_0\rangle$ is a state with norm 1)
More succinctly, if there exists a unitary operator $U$ such that $[H,U]=0$ and $U|\phi_0\rangle\neq e^{i\theta}|\phi_0\rangle$ for any phase $\theta$ then we have ground state degeneracy.
In the example you have given we see that the matrix elements in the basis given $\{|a\rangle,|b\rangle,|c\rangle\}$ is $$H = \begin{pmatrix}1&0&0\\0&1&0\\0&0&2\end{pmatrix}$$ from which we see there exists a unitary operator, with matrix elements $$U = \begin{pmatrix}0&1&0\\1&0&0\\0&0&1\end{pmatrix}$$ which commutes with $H$ and is non-trivial in the ground-state space.
Proof that non-existence of $U$ implies non-degenerate ground state:
Assume $\nexists U$ s.t. $\{[H,U]=0 ~~\mbox{and}~~ U|\phi_0\rangle \neq e^{i\theta}|\phi_0\rangle\}$
Now, for every state $|a\rangle$ and $|b\rangle$, $\exists U_{ab}$ which is unitary that takes us from $|a\rangle\rightarrow|b\rangle$. We are interested in the operator that take us from $|\phi_0\rangle$ to any $|a\rangle$ in our Hilbert space (which obviously includes all possible ground states), which we denote by $U_{a0}$. This means that any state $|a\rangle$ can be written as $|a\rangle = U_{a0}|\phi_0\rangle$. By our starting assumption $U_{a0}$ either satisfies
$$(1)~~~~~~ [H,U_{a0}]\neq 0,~~~~~~~~\mbox{or}~~~~~~~~(2)~~~~~U_{a0}|\phi_0\rangle = e^{i\theta}|\phi_0\rangle$$
If (1), then we have $$H|a\rangle = H U_{a0}|\phi_0\rangle \neq U_{a0}H|\phi_0\rangle = E_0|a\rangle~~~\implies~~~H|a\rangle \neq E_0|a\rangle$$
and so $|a\rangle \neq |\phi_0\rangle$ is not a ground state.
If (2), then $|a\rangle = e^{i\theta}|\phi_0\rangle$ and so $|a\rangle$ and $|\phi_0\rangle$ represent the same state.
Thus the non-existence of $U$ implies the non-existence of a second ground state and thus non-degeneracy.
Best Answer
Meng Cheng pretty much answered the question already, but it might still be helpful to elaborate on the meaning of "too hard to calculate $e^{-\tau H}$ for a finite $\tau$".
What hardness refers to is the computational complexity of the problem/proposed numerical algorithm. That is to say, given a problem of size $N$ (e.g. here in the paper referenced $N$ is the number of spins wished to simulate), how much memory and how much time do you need.
Why this is important is because we typically want to solve problems in which $N$ is big: obviously, computing $e^{-\tau H}$ for say a 1024x1024 Hamiltonian $H$ describing 10 spins is not a problem -- just type exp(H) in MATLAB or Mathematica. But ask MATLAB to exponentiate $H$ for 100 spins and your computer will crash.
Since exponentiation of a Hamiltonian typically results in a dense matrix, this has bad memory scaling: exponential in $N$. For scaling in time, see: https://mathoverflow.net/questions/239073/what-is-the-time-complexity-of-the-matrix-exponential
Therefore, naively trying to compute $e^{-\tau H}$ is folly.
The key insights that the paper uses are:
Exponential over finite time can be written as product of many exponentials over small times $e^{-\tau H} = e^{-d \tau H}e^{-d \tau H}\cdots e^{-d \tau H}$ (this is not an approximation).
For small time exponentials $e^{-d \tau H}$, it can be approximated (Trotter decomposition) as in Eq. 4 of the paper: $e^{-d\tau H} = e^{-d\tau H_z} e^{-d\tau H_y} e^{-d\tau H_x} + O(d\tau^2)$
Naively you might think, $H_z, H_y, H_x$ are still big matrices, so exponentiating them is still hard, like in the general case. But what's nice about each $H_z, H_y, H_x$ is that they are each made of sums of local terms which commute. So the full matrix $e^{-d\tau H_z}$ is itself simply made of products of local exponentiations $e^{-d\tau H_{z,local}}$ where $H_{z,local}$ only acts on 4 sites (assuming a nearest-neighbor model) and is thus easy to do.
Furthermore, note that we intend to apply $e^{-\tau H}$ to a wavefunction. So we don't actually have to form the full $e^{-d \tau H_z}$ etc. which would be humongous. Instead we just apply $e^{-d\tau H_{z,local}}$ to each local part of the wavefunction. This is very fast and memory efficient!