This Asymptote code seems to work perfectly and for any $t$ in your range the estimate uses $Cd$ operations and is a guaranteed upper bound though I am not sure whether $C$ is small enough for you (I hope it is).
If you try to run it online on http://asymptote.ualberta.ca/ , the final pause() command should be removed. The output gives the value of $t$ (I run it over powers of 2 to make the comparison quick), the truth and the quick upper bound after it. All explanations tomorrow.
int d=300;
srand(seconds());
real p=0.6+2*unitrand();
real[] h,A;
real[][] C;
for(int k=0;k<d;++k) h[k]=exp(-p*log(k+1))*(1+2*unitrand());
h=reverse(sort(h));
real a=1/(2*max(h)+sum(h));
for(int k=0;k<d;++k)
{
C[k]=new real[];
for(int l=0;l<d;++l) {C[k][l]=a^2*h[k]*h[l];}
real u=a*h[k];
A[k]=2*(u-u^2); C[k][k]+=1-2*(u-u^2);
}
real S(real[][] A)
{
real s=0; for(int k=0;k<d;++k) for(int l=0;l<d;++l) s+=A[k][l];
return s;
}
real[][] U=copy(C);
for(int kk=0;kk<19;++kk)
{
int t=2^kk;
write(t);
real la=A[d-1]-a^2*h[d-1]^2;
for(int p=0;p<4;++p)
{
real s=-1,ss=0;
for(int k=0;k<d;++k) {s+=a^2*h[k]^2/(A[k]-la); ss+=a^2*h[k]^2/(A[k]-la)^2;}
la-=s/ss;
}
pair r1=(-3.45888846213541,1.77794875224373), r2=(-0.541111537864589,5.00955180872487),
d1=(0.0960769769474388,0.281952388994648), d2=(-0.762743643614106,-0.350175117369014);
write(S(U));
real s=0; pair ss=(0,0),sss=(0,0);
pair r=-1/r1;
for(int k=0;k<d;++k) {ss+=a*h[k]/(1+t*r*(A[k]-la)); sss+=(a*h[k])^2/(1+t*r*(A[k]-la));}
pair al=ss/(1-t*r*sss);
pair[] Q;
for(int k=0;k<d;++k) {Q[k]=-(2/d1/r1)*((1+t*r*al*a*h[k])/(1+t*r*(A[k]-la)));}
r=-1/r2;
ss=(0,0); sss=(0,0);
for(int k=0;k<d;++k) {ss+=a*h[k]/(1+t*r*(A[k]-la)); sss+=(a*h[k])^2/(1+t*r*(A[k]-la));}
pair al=ss/(1-t*r*sss);
for(int k=0;k<d;++k) {Q[k]-=(2/d2/r2)*((1+t*r*al*a*h[k])/(1+t*r*(A[k]-la)));}
for(int k=0;k<d;++k) s+=Q[k].x^2;
s*=(1-la)^t;
write(s);
write("************");
U*=U;
}
pause();
Now the explanations.
To the best of my understanding, you are aiming at finding $\langle(1-V)^t w,w\rangle$ with low relative error where $w$ is the vector of all $1$'s and $V=A-h\otimes h$ where $A=\text{diag}(a_j)$ and $h=(h_j)$ with $a_j,h_j>0$ ($j=1,\dots,d$).
The first observation is that we can factor out $(1-\lambda)$ where $\lambda$ is the least eigenvalue of $V$. It is the $\lambda$ that makes possible to have an equality in the positive definiteness Cauchy-Schwarz inequality $\sum_j{a_j-\lambda}x_j^2\ge\left(\sum_j h_jx_j\right)^2$ which is the least root of the equation
$$
\sum_j\frac{h_j^2}{a_j-\lambda}=1\,.
$$
I'm finding it by Newton iterations (I used 4 in the above program but if you want to go really high in $t$, 6 would be better). Note that it is crucial to get the correct root, so one should start with the initial approximation $a_{j_0}-h_{j_0}^2$ where $a_{j_0}=\min_j a_j$. It is tempting to start with $0$, but that will often throw you in between $a_j$ and result in a wrong root, which will ruin everything.
Now we have $V=\lambda I+V'$ where $V'=A'-h\otimes h$, $A'=\text{diag\,}(a_j-\lambda)$ is a SPD matrix with the least eigenvalue $0$ and the corresponding unit eigenvector $e_0$ that has positive entries (proportional to $\frac{h_j}{a_j-\lambda}$). This
$$
(I-V)^t=(1-\lambda)^t(1-\tfrac 1{1-\lambda}V')^t\approx (1-\lambda)^t\exp(-\frac t{1-\lambda}V')
$$
and $\approx$ is an upper bound.
So, from now on, I'll remove the primes, denote $a_j-\lambda$ by $a_j$, $V'$ by $V$, and $\frac{t}{1-\lambda}$ by $t$ (I didn't divide $t$ by $(1-\lambda)$ in the code above because $\lambda$ was very small in your setting and you do not need it when finding $e^{-tV}$ but if you are after the discrete case, it helps a bit).
Theoretically, $\langle e^{-tV}w,w\rangle$ (BTW, I recommend computing it after this renormalization because the geometric progression $(1-\lambda)^t$ doesn't affect the relative error but can easily result in mantissa overflow) is just $\langle w,e_0\rangle^2+\sum_{j\ge 1}e^{-\gamma_j t}\langle w,e_j\rangle^2$ where $\gamma_j$ and $e_j$ are the remaining eigenvalues and eigenvectors. The idea is to replace it with $\langle F(tV)w,w\rangle$ where $F(t)\ge e^{-t}$ is a good approximation to $e^{-t}$ (we will discuss below how good we want it) and $\langle F(tV) w, w\rangle$ is reasonably quick to compute.
All I really know how to compute quickly is $w(z)=(I+zV)^{-1}w$. This calls for solving the system of equations
$$
w(z)_i (1+za_i)-z\langle w(z),h\rangle h_i=1\,.
$$
Let $\alpha=\langle w(z),h\rangle$. Then $w(z)_i=\frac{1+z\alpha h_i}{1+z a_i}$ and $\alpha$ can be found from
$$
\sum_i \frac{1+z\alpha h_i}{1+z a_i}h_i=\alpha\,,
$$
i.e., $\alpha=\frac{s_1}{1-zs_2}$ where $s_q=\sum_i \frac{h_i^q}{1+z a_i}$, $q=1,2$.
Thus we can afford rational functions of not too high degree. This is the key. We can now use $F(t)=P(t/2)^{-2}$ where $P$ is the partial sum of the Taylor series for $e^t$ (I used $P(t)=1+t+t^2/2+t^3/6+t^4/24$ in the above program because I was aiming at the range of hundreds (Asymptote will just run out of memory if you try to make a $d\times d$ matrix with $d=10^4$ and it is also about 100 times slower than C++, because it is an interpreter, not compiler; for $d\asymp 10^4$ I would rather recommend degree 8).
If $r_m$ are the (complex!) roots of $P$ (which have to be found once with high precision; I actually found the roots of $P(t/2)$ in the upper half-plane, those are the mysterious $r1$ and $r_2$ in the code) and $d_m=P'(r_m)$, then we have the representation
$$
\frac 1{P(t)}=\sum_m\frac {d_m}{t-r_m}=-\sum_m \frac 1{r_md_m}\frac 1{1+(-1/r_m)t}=-2\Re {\sum_m}' \frac 1{r_md_m}\frac 1{1+(-1/r_m)t}
$$
where ${\sum}'$ is taken over the roots in the upper half-plane (so you have to run the above computation just for 2 roots for degree 4, 4 roots for degree 8, etc.).
Thus we get
$$
P(tV/2)^{-1}w=-2\sum'_m \Re[(r_md_m)^{-1}w(-r_m^{-1})t/2]
$$
which is readily computable as above, and
$\langle P(tV/2)^{-2}w,w\rangle=\|P(tV/2)^{-1}w\|^2$ (one more summation of squares of entries).
Now it is time to discuss precision.
Since $\sum_{j\ge 0}\langle w,e_j\rangle^2=\|w\|^2=d$, we have
$$
\langle e^{-tV}w,w\rangle=\sum_{j\ge 1}\langle w,e_j\rangle^2\left[\tfrac {\langle w,e_0\rangle^2}{d-\langle w,e_0\rangle^2}+e^{-t\gamma_j}\right]
$$
and
$$
\langle F(t)w,w\rangle-\langle e^{-tV}w,w\rangle=
\sum_{j\ge 1}(F(t\gamma_j)-e^{-t\gamma_j})\langle w,e_j\rangle^2\,.
$$
Note that all terms are positive, so the relative error in the sum is at most the maximum of the relative errors in the individual terms. Note also that, since $e_0$ has positive entries, we have $\langle w,e_0\rangle\ge 1$. Thus we are interested in
$$
E=\max_{t>0}\frac{F(t)-e^{-t}}{\frac 1{d-1}+e^{-t}}=
max_{t>0}\frac{P(t/2)^{-2}-e^{-t}}{\frac 1{d-1}+e^{-t}}\,.
$$
This can be investigated theoretically, but I'll just make a small table (the first number is the degree of the polynomial and the second is $d$) so that you can see that for $n=4, d=100$ I have only 10% error in the worst case scenario (your data can never go that high) while for $d=10^4$, with $n=4$ I can theoretically go 150% above the truth while $n=8$ gives you guaranteed 5% precision at the expense of 2 extra roots. If you can afford $n=14$, you can go to ten millions range and declare your computations "precise".
Edit: The $\lambda$-story.
The estimate $E=E(d,n)$ for the relative error above is uniform in all $t>0$, but it is based of the assumption that the lowest eigenvalue $\gamma_0$ is exactly $0$, so not only the term $\langle w,e_0\rangle^2e^{-\gamma_0 t}$ is evaluated exactly by $\langle w,e_0\rangle^2 F(t)$, but it also stays fixed independently of $t>0$. The relative error bound is based not on the idea that $F(t)$ approximates $e^{-t}$ with small relative error on the whole semi-axis (that is impossible and you have a blow-up beyond $t=40$ or so for $n=16$) but on the idea that both $F(t)$ and $e^{-t}$ become negligibly small compared to $\frac 1{d-1}$ when the relative error between them becomes large.
In actual computation, you'll not be able to keep $\gamma_0$ exactly at $0$ (the computation of $\lambda$ has some error), so it is important to understand what happens if it is small. The answer is that as long as $|\gamma_0|t<0.01$, say, nothing changes in the bounds; just $1$ gets replaced by $0.99$ in the formula for $E$. However, if this product is large, the whole thing falls apart dramatically because you accumulate a huge relative error in the term that was presumed to be evaluated exactly and the numerator in $\frac 1{d-1}$ will also become $e^{-\gamma_0 t}$. So, you are theoretically guaranteed to have the relative error as claimed only if you use the $\lambda I$ reduction and compute $\lambda$ with the error at most $0.01 t^{-1}$. Otherwise you'll be outputting junk for large $t$. On the other hand, the theoretical guarantee (once you just check that its assumptions hold) is a sure one, so you do not need to check that its conclusion holds directly.
Moral: make sure that the error in $\lambda$ is much smaller than $1/t$ and you will not need to check anything else after that, provided that there is no stupid error in the code, which would reveal itself on small data.
2 10 0.17245166347402
2 100 0.961455640371778
2 1000 4.28015125191342
2 10000 19.7816887991447
2 100000 101.301563701775
2 1000000 570.64054529331
2 10000000 3464.50750307103
*********
4 10 0.0162756483432251
4 100 0.106015462775259
4 1000 0.429607071204345
4 10000 1.4053481903916
4 100000 4.41442731700713
4 1000000 14.7171384159884
4 10000000 54.2101332254995
*********
6 10 0.00166712645387837
6 100 0.0139026251293181
6 1000 0.0705786452201884
6 10000 0.2484898912299
6 100000 0.708166083384971
6 1000000 1.85420759122078
6 10000000 4.87739303482548
*********
8 10 0.00016984164765101
8 100 0.00166787956056493
8 1000 0.0110935381765316
8 10000 0.0492442082940087
8 100000 0.159219723419999
8 1000000 0.420853472933991
8 10000000 1.00146598253247
*********
10 10 1.73531792069359e-05
10 100 0.000183806352715792
10 1000 0.00151002712819166
10 10000 0.00873278982007405
10 100000 0.0352404205756958
10 1000000 0.107570103457829
10 10000000 0.271655462239004
*********
12 10 1.78880885047298e-06
12 100 1.94620128830683e-05
12 1000 0.000180991851874676
12 10000 0.00131944915308362
12 100000 0.00681706198345839
12 1000000 0.0256053500194781
12 10000000 0.0750096474609068
*********
14 10 1.86072838564969e-07
14 100 2.04069133751507e-06
14 1000 2.00507453637578e-05
14 10000 0.000172432808234548
14 100000 0.00112204219089126
14 1000000 0.0052891118636613
14 10000000 0.0187925172629428
