Return to Question

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

I'm looking to do this in $O(d)$ time by specializing to the case of convergent matrices. There's a suggestion that this is doable in Section 3 of this paper so I'm looking for hints on the best way to implement this (or a cheap alternative) in Mathematica.

Example below uses generating function approach to evaluate powers of convergent diagonal matrices. It uses approximation below, then calls LaplaceTransform on integral obtained after IntegrateChangeVariables

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

It should print

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

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];
   
   Assert[Head[int] == Inactive[Integrate]];
   {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]];

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

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];
   
   Assert[Head[int] == Inactive[Integrate]];
   {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[5]

I'm looking to compute the following quantity for different values of $k$, 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 to do this in $O(d)$ time by specializing to the case of convergent matrices. There's a suggestion that this is doable in Section 3 of this paper so I'm looking for hints on the best way to implement this (or a cheap alternative) in Mathematica.

Example below uses generating function approach to evaluate powers of convergent diagonal matrices. It uses approximation below, then calls LaplaceTransform on integral obtained after IntegrateChangeVariables

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

It should print

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

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];
   
   Assert[Head[int] == Inactive[Integrate]];
   {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]];

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^d)$$

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 to do this in $O(d)$ time by specializing to the case of convergent matrices. There's a suggestion that this is doable in Section 3 of this paper so I'm looking for hints on the best way to implement this (or a cheap alternative) in Mathematica.

Example below uses generating function approach to evaluate powers of convergent diagonal matrices. It uses approximation below, then calls LaplaceTransform on integral obtained after IntegrateChangeVariables

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

It should print

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

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];
   
   Assert[Head[int] == Inactive[Integrate]];
   {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]];

I'm looking to compute the following quantity for different values of $k$, 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 to do this in $O(d)$ time by specializing to the case of convergent matrices. There's a suggestion that this is doable in Section 3 of this paper so I'm looking for hints on the best way to implement this in Mathematica.

Example below uses generating function approach to evaluate powers of convergent diagonal matrices. It uses approximation below, then calls LaplaceTransform on integral obtained after IntegrateChangeVariables

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

It should print

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

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];
   
   Assert[Head[int] == Inactive[Integrate]];
   {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]];

I'm looking to compute the following quantity for different values of $k$, 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 to do this in $O(d)$ time by specializing to the case of convergent matrices. There's a suggestion that this is doable in Section 3 of this paper so I'm looking for hints on the best way to implement this (or a cheap alternative) in Mathematica.

Example below uses generating function approach to evaluate powers of convergent diagonal matrices. It uses approximation below, then calls LaplaceTransform on integral obtained after IntegrateChangeVariables

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

It should print

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

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];
   
   Assert[Head[int] == Inactive[Integrate]];
   {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]];