I have non-negative $d\times d$ matrices $A$, $B$ and need a tractable way to compute the sum of all entries of $\exp(-t(A-B))$ where $A$ is diagonal and $B$ symmetric rank-$1$. IE


$$f(t)=\langle\exp(-t(A-B))\rangle$$


Where $t>0$ and $\langle M\rangle$ represents sum of all entries of $M$.


**What expansion would be useful to estimate $f(t)$?**


I know that norms of $A$, $B$ and $A-B$ are less than one, $A-B$ is positive semi-definite and norm of $B$ is relatively small compared to norm of $A$.


I can only afford to perform $c\times d$ operations for small values of $c$.
Golden-Thompson [inequality](https://terrytao.wordpress.com/2010/07/15/the-golden-thompson-inequality/) gives efficient to compute and tight bound for $\operatorname{Tr}$, but similar factoring is loose for $\langle\cdot\rangle$. Computing $f(t)$ would extend Velikanov's [approach](https://arxiv.org/abs/2206.11124) for computing SGD error after $t$ steps to small values of $t$.


[Notebook](https://www.wolframcloud.com/obj/yaroslavvb/nn-linear/forum-bounding-rank1.nb)



Approximating sum of entries of $\exp(A-B)$ for diagonal $A$ and rank-$1$ $B$?
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
 

```

