Speed up minimization of a high-degree polynomial

Question

For fixed $s$ and $n$, I need to find $\alpha$ which minimizes the following expression

$$h (I-2\alpha H + 2\alpha^2 H + \alpha^2 \mathcal{H})^s x0$$ where

$$ \begin{array}{lll} h&=&(1,\frac{1}{2},\frac{1}{3},\ldots,\frac{1}{n}) \\ H&=&\text{diag}(h)\\ \mathcal{H}&=& h\otimes h\\ x0&=& (1,1,\ldots,1)^T \end{array} $$

I got the straightforward solution by using NMinimize, but it already takes 7 seconds for $n=5$, $s=10$, I'm wondering if there are Mathematica tricks to make this feasible for larger $n$ and $s$

CircleTimes = KroneckerProduct;
bestRate[n_, numSteps_] := Module[{},
   ii = IdentityMatrix[n] // N;
   h = Table[i^-1, {i, 1, n}] // N;
   hh = h\[CircleTimes]h;
   H = DiagonalMatrix[h];
   M = ii - 2 alpha H + 2 alpha^2 H.H + alpha^2 hh;
   x0 = ConstantArray[1., {n}];
   xn = MatrixPower[M, numSteps].x0;
   alpha /. Last@NMinimize[xn.h, alpha]
   ];
With[{n = 2},
 Table[bestRate[n, s], {s, 1, 20}]
 ]

PS: this is the task of finding optimal rate for Least mean squares filter which runs for $s$ steps, with observations coming from Normal(0,H)

Answer 1

Exact root-finding is very fast here:

bestRate[n_, numSteps_] := Module[{ii,h,hh,H,M,x0,xn,α},
  ii = IdentityMatrix[n];
  h = 1/Range[n];
  hh = h~KroneckerProduct~h;
  H = DiagonalMatrix[h];
  M = ii - 2 α H + 2 α^2 H . H + α^2 hh;
  x0 = ConstantArray[1, {n}];
  xn = MatrixPower[M, numSteps] . x0;
  α /. Last@Solve[Expand[D[xn . h, α]] == 0, α, Reals]]

With[{n = 2},
  Table[bestRate[n, s], {s, 1, 20}]]

(*    {10/33, 0.30446, 0.30617, 0.308113, 0.310232, 0.312465,
       0.314754, 0.317042, 0.319286, 0.321449, 0.323507, 0.325446,
       0.327256, 0.328938, 0.330493, 0.331927, 0.333246, 0.33446,
       0.335576, 0.336602}                                            *)

The results are given as Root objects. Here's an instructive answer on how to work with these.

With[{n = 2},
  DiscretePlot[bestRate[n, s], {s, 1, 100}]]

I think the asymptotic $n=2$ answer for $s\to\infty$ is

Root[32 - 288 # + 974 #^2 - 1449 #^3 + 776 #^4 &, 1]
(*    0.356182    *)

Answer 2

For me, using Simplify when calculating xn greatly reduced the time for minimization. It was clear from using AbsoluteTiming that the majority of time was spent during minimization for your use case.

Original code ran in 6.4 seconds with n=5 and s=10. For the same parameters, the following code ran in 0.84 seconds for me. I'm sure there are more quick improvements that are possible by just simplifying the function to optimize.

bestRate1[n_, numSteps_] := Module[{}, ii = IdentityMatrix[n] // N;
   h = Table[i^-1, {i, 1, n}] // N;
   hh = h\[CircleTimes]h;
   H = DiagonalMatrix[h];
   M = ii - 2 alpha H + 2 alpha^2 H.H + alpha^2 hh;
   x0 = ConstantArray[1., {n}];
   xn = Simplify[MatrixPower[M, numSteps].x0];
   alpha /. Last@NMinimize[xn.h, alpha]];

Things that were tried but didn't do much

1. Using ParallelTable. I doubt this is useful for you since you probably don't plan on setting n=1000.
2. Using FullSimplify at various stages (i.e. M, xn). Took too long.
3. Using Refine at various stages (i.e. M, xn). Probably didn't do any worthwhile simplification from the perspective of the optimizer.

Another thing you could try out is using a different minimization method. NelderMead and DifferentialEvolution seem to give similar results. But for the values you had, SimulatedAnnealing and RandomSearch showed marginal improvement.

Answer 3

I would black-box the objective function so that the matrix powering does not happen symbolically. Also I would use a vector as third argument of the matrix power step. These steps really do not help much for small s and n but can make a big difference as one or both grow.

makeMat[n_, alf_] := Module[
  {h, H},
  h = 1./Range[n];
  H = DiagonalMatrix[h];
  IdentityMatrix[n] - 2 alf H + 2 alf^2 H . H + 
   alf^2*KroneckerProduct[h, h]
  ]

objfunc[mat_?(MatrixQ[#, NumberQ] &), s_Integer] := 
 With[{n = Length[mat]}, 
  MatrixPower[mat, s, ConstantArray[1., {n}]] . (1/Range[n])]

bestRate2[n_, s_] := Module[
  {mat, alf},
  mat = makeMat[n, alf];
  alf /. Last[NMinimize[objfunc[mat, s], alf]]
  ]

Test at larger values of s:

In[39]:= AbsoluteTiming[
 With[{n = 2}, Table[bestRate2[n, s], {s, 26, 50}]]]

(* Out[39]= {1.79803, {0.341282, 0.341872, 0.34242, 0.342931, 0.343407, 
  0.343852, 0.344269, 0.344659, 0.345026, 0.34537, 0.345695, 0.346002,
   0.346291, 0.346565, 0.346824, 0.34707, 0.347303, 0.347525, 0.34774,
   0.347937, 0.348129, 0.348311, 0.348488, 0.348655, 0.348811}} *)

The original version took 33.3 seconds on the same machine.