In[]:=
(*deployswithcanonicalname*)deploy
Thu 16 Mar 2023 19:23:41
Upper bounds on exp(A-B)
Upper bounds on exp(A-B)
Numerically check upper bounds on realistic example.
Mathoverflow question: https://mathoverflow.net/questions/442888/approximating-sum-of-entries-of-expa-b-for-diagonal-a-and-rank-1-b
Mathoverflow question: https://mathoverflow.net/questions/442888/approximating-sum-of-entries-of-expa-b-for-diagonal-a-and-rank-1-b
In[]:=
outer2[h_]:={h}.{h};(*Returnsminibatchoperatoras{diag,rank1}*)minibatchOperatorFactored[h_,α_?NumericQ,b_]:=Module{d=Length[h],H,Hm,ones,ii},ones=ConstantArray[1.,d];Ifb==∞,{ones-2αh+(h*h),0*ones},ones-2αh+(b+1)(h*h),αh(*sameas+(h*h),αh*);minibatchOperator[h_,α_?NumericQ,b_]:=Module[{d,z},{d,z}=minibatchOperatorFactored[h,α,b];DiagonalMatrix[d]+outer2[z]];stepL1[h_,b_]:=With{r=Total[h]/Max[h]},Ifb==∞,2Max[h],+1;specRadius[mat_]:=Max@Eigenvalues[mat,1];stepL2[h_,b_]:=With{hh=N@Outer[Times,h,h],H=N@DiagonalMatrix[h]},Ifb==∞,2Norm[H],;(*Setuprealisticproblem.a+b/.subreplacesA+Bwithrealisticdiag+rank1*)traceNormalize[h_]:=h/Total[h];d=100;xvals=Range[1,d];ones=ConstantArray[1.,d];p=1;g[x_]=;h=g/@xvals;h=h/Total[h];b=1;(*batchsize*)s=d;(*stepcount*)α=0.7*stepL2[h,b];α=0.5*stepL2[h,b];α=0.9999*stepL2[h,b];α=N[1/Tr[h]];(*stepsize*)ii=IdentityMatrix[d];A=DiagonalMatrix2αh--;(*diagonalterm*)B=-{h}.{h};(*rank1term*)On[Assert];T=minibatchOperator[h,α,b];Assert[ii-(A+B)==T];true=Tr@MatrixPower[ii-(A+B),s];eval[q_,label_]:={q-true,Log[Abs[q-true]],q,label};vals={eval[Tr@MatrixPower[ii-(A+B),s],"exact"],eval[Total[],"gd"],eval[Total[],"diagonal"],eval[Total[Exp[-sDiagonal[(A+B)]]],"diagonal exp"],eval[Tr[MatrixExp[-s(A+B)]],"full exp"],eval[Tr[MatrixExp[-s(A)].MatrixExp[-s(B)]],"Golden-Thompson"]};SF=StringForm;Print[SF["Traces for d=``, p=``, b=``, s=``",d,p,b,s]];TableForm[vals[[All,;;-2]],TableHeadings->{vals[[All,-1]],{"error","logerror","val"}}]Print[SF["Errors for d=``, p=``, b=``, s=``",d,p,b,s]];(*testcasefordiagonalcovariance,*)On[Assert];diag[mat_]:=DiagonalMatrix[Diagonal@mat];Assert[b==1];step[c_]:=(+)c+hTotal[hc];true=Total[Nest[step,ones,s]];(*Evaluatesumofallentries*)err[mat_]:=Total[mat.ones];eval[mat_,label_]:=With[{q=err[mat]},{q-true,Log[Abs[q-true]],q,label}];(*trunatedbinomialexpansion*)Clear[a,b];term={b}~Join~Table[a,s-1];terms=Table[RotateRight[term,ss],{ss,0,s}];evalTerm[term_]:=Dot@@term/.{a->ii-A,b->-B}order0=evalTerm[Table[a,s]];order1=Total[evalTerm/@terms];order1mat=order0+order1;order1approx=order0+evalTerm[terms[[s/2]]]*s;order2approx=order1approx+evalTerm[RandomSample[{b,b}~Join~Table[a,s-2],s]];order3approx=order2approx+evalTerm[RandomSample[{b,b,b}~Join~Table[a,s-3],s]];
2
α
2
α
b
1
b
2
(1-αh)
2
α
1
b
1
b
2
Total[h]
b
(b+1)
r
2b
specRadius[(b+1)H+hh.Inverse[H]]
-p
x
2
α
2
h
2
α
b
2
h
2
α
b
2s
(1-αh)
s
Diagonal[ii-(A+B)]
2s1
(1-h)
2
(1-αh)
2
α
2
h
2
α
s!
2!(s-2)!
s!
3!(s-3)!
Traces for d=100, p=1, b=1, s=100
Out[]//TableForm=
error | logerror | val | |
exact | 0. | Indeterminate | 40.4378 |
gd | -0.193834 | -1.64076 | 40.244 |
diagonal | -0.0606168 | -2.80318 | 40.3772 |
diagonal exp | 0.0711606 | -2.64282 | 40.509 |
full exp | 0.131662 | -2.02751 | 40.5695 |
Golden-Thompson | 4.76912 | 1.56216 | 45.2069 |
Errors for d=100, p=1, b=1, s=100
Check constraints of “Approximating sum” question
Check constraints of “Approximating sum” question
Must run previous cell.
https://mathoverflow.net/questions/442888/approximating-sum-of-entries-of-expa-b-for-diagonal-a-and-rank-1-b
https://mathoverflow.net/questions/442888/approximating-sum-of-entries-of-expa-b-for-diagonal-a-and-rank-1-b
In[]:=
Bnn=-B;(*non-negativerank-1matrix*)On[Assert];Assert[And@@Thread[{Norm[A],Norm[Bnn],Norm[A-Bnn]}<1]]Assert[And@@(PositiveSemidefiniteMatrixQ/@{A,Bnn,A-Bnn})]nnMatrixQ[mat_]:=And@@Thread[Flatten[mat]>=0]Assert[nnMatrixQ[A]]Assert[nnMatrixQ[Bnn]]vals={eval[Nest[T.#&,T,s-1],"exact"],eval[,"gd"],eval[MatrixPower[diag[ii-(A+B)],s],"diagonal"],eval[MatrixPower[ii-A,s],"order0"],eval[order1approx,"order1-approx"],eval[order2approx,"order2-approx"],eval[order3approx,"order3-approx"],eval[order1mat,"order1"],eval[MatrixExp[-s(A+B)],"full exp"],eval[MatrixExp[-s(A-Bnn)],"full exp2"],eval[MatrixExp[-s(A)],"bound"],eval[MatrixExp[-sdiag[(A+B)]],"diag exp"],eval[MatrixExp[-s(A)].MatrixExp[-s(B)],"Golden-Thompson"]};TableForm[vals[[All,;;-2]],TableHeadings->{vals[[All,-1]],{"error","logerror","val"}}]terms=Table[RotateRight[term,ss],{ss,0,s}];ListPlot[Total[evalTerm[#],2]&/@terms,PlotLabel->"values of ||AAABAAA|| for B in k'th position",AxesLabel->{"k","value"}]
2s
(1-αh)
Upper bounding diagonal + rank1
Upper bounding diagonal + rank1