# Efficient powers of DPR1 matrices

I'm looking to compute the following quantity for $$A$$ where $$A$$ is a convergent positive definite $$d\times d$$ diagonal + rank1 matrix (DPR1).

$$f(s)=\operatorname{Tr}(A^s)$$

Earlier answer answer by Henrik Schumacher shows how to do this in Mathematica in $$O(d^2)$$ time for an arbitrary positive definite DPR1 matrix, and provides a nice package for this.

I'm looking for the most elegant way to do this in Mathematica in $$O(d)$$ time by specializing to the case of convergent matrices. The math for this technique is described here

Example below uses LaplaceTransform to find $$\operatorname{Tr}(A^s)$$ when $$A$$ is diagonal. It relies on following approximation.

$$\sum_i (1-h_i)^s\approx \int_i \exp(-sh_i)$$

d = 5;
h[i_] = i^-1;
hvals = Table[h[i], {i, 1, d}];
diag[s_] := Tr[MatrixPower[DiagonalMatrix[1 - hvals], s]]/d;
diagInt = Inactive[Integrate][Exp[-s h[i]]/d, {i, 1, d}];

(* changes integration limits (a,b) -> (0,\[Infinity]) *)
expandIntegrationLimits[int_] :=
Module[{body, limitSpec, var, limits, newLimitSpec, interval},
(* returns expression that is zero when var is outside of (min,
var) *)
interval[var_, {min_, max_}] =
HeavisideTheta[var - min] - HeavisideTheta[var - max];

{body, limitSpec} = List @@ int;
{var, limits} = {First[limitSpec], Rest[limitSpec]};
newLimitSpec = {var, 0, \[Infinity]};
Inactive[Integrate][body*interval[var, limits], newLimitSpec]
];

(* turn into Laplace transform form *)
int = IntegrateChangeVariables[diagInt, hi, hi == h[i]];
int = expandIntegrationLimits[int];
diagApprox[s_] = -LaplaceTransform[First@int/Exp[-hi s], hi, s];

s0 = 5;
Print["Exact value diag: ", diag[s0] // N]
Print["Integral approximation: ",
Activate[diagInt /. {Integrate -> NIntegrate, s -> s0}]]

Print["symbolic Laplace: ", N[diagApprox[s0]]];

(* normalize h to ensure DPR1 matrix is convergent *)
hvals = hvals/Total[hvals];
dpr[s_] :=
Tr[MatrixPower[
DiagonalMatrix[1 - hvals] + {hvals}\[Transpose] . {hvals}, s]]/
d;
s0 = 5;
Print["Exact value dpr1: ", dpr[s0] // N];
Print["Approx value dpr1: ", dprApprox[s0]];


should print

Exact value diag: 0.145584
Integral approximation: 0.148296
symbolic Laplace: 0.148296
Exact value dpr1: 0.528644
Approx value dpr1: dprApprox