How to make the computation of a minimum efficient?

Question

I have a $80\times 80$ matrix $A$ and a vector $P$ and I want to compute efficiently the following minimum $$ \min\{n\in \Bbb N :\|A^nP\|_\infty < 1\} $$ Then I defined my matrix and vector in mathematica with the following code

A = SparseArray[{{1, 30} -> 1/2, {i_, j_} /; i - j == 1 -> 1}, {80, 
80}]; 
P = Table[8 10^9 2^(-80) Binomial[80, k], {k, 0, 79}]

and to compute the above minimum I tried the following code

Minimize[{
    Norm[MatrixPower[A, n].P, Infinity] < 1, n ∈ Integers, n > 0
  }, n]

I tried different versions of the above code using NMinimize instead of Minimize, or N@Norm or Max instead of just Norm, but every time my computer stay calculating for more than five minutes without giving a result.

However using directly Norm and MatrixPower and different values for n I get manually the result for $n=801$ in less than 30 seconds, so I assume there is a way to compute the above efficiently, and this is the reason why I had opened this question.

Someone knows some way to compute the above fast?

---

EDIT: the original question is above. I had wrote wrongly the code to find the minimum value of $n$ such that $\|A^nP\|_\infty <1$. However I dont know how to fix it. Can someone help me?

I tried the code

Minimize[{n, 0 < Max[MatrixPower[A, n].P] <= 1 && n > 0}, 
n ∈ Integers]

but it still stay in the computation for a long time, so I cancelled it.

Answer 1

Never underestimate the power and ease of whipping together a binary search.

Your function is monotonically non-increasing, as you noted. So you can use a binary search to robustly find the value you want.

binSearch[left_, right_] := Module[{l = left, r = right, middle = Floor[(right + left)/2]},
    While[(r - l > 1), 
          (
           If[Max[MatrixPower[A, middle].P] - 1 < 0, r = middle, l = middle];
           middle = Floor[(r + l)/2]; 
           (* watch the convergence *)
           Print["{" <> ToString@l <> "," <> ToString@r <> "}"];
           )
       r]]

Pick values far enough apart to ensure the answer lies between them

binSearch[1, 10000]
(* 801 *)

This will converge in $\log_{2}10000 = 13.2877\approx 13$ iterations. Not too shabby.

EDIT

As an OBTW, it's interesting to take a look at the form of the matrix $A^{n}$

Manipulate[MatrixPlot[Log@MatrixPower[A, k], Mesh -> All], {k, 1, 1000, 1}]

From eyeballing it, you can show that the diagonals are just powers $2^{-k}$ for integers $k$, that you can replace P with

Psimp = {0, 0, 0,...,P[[29]], P[[30]], 0, 0,..., 0},

that the function is monotonic nonincreasing, and that after a transient phase as you move through n,all of the unique values in your function can be duplicated simply...

brute= Table[Max[MatrixPower[A, k].P], {k, 80, 201, 1}] // Intersection // Reverse;

elegant = Flatten@Table[{P[[30]]/2^(k - 1), P[[29]]/2^(k - 1)}, {k, 2, 6}];

brute == elegant
(* true *)

With a little bit of fiddling, a simple closed form expression is within reach.

Answer 2

Your function is non-trivial. Here is a log-plot for integer values of $n$:

DiscretePlot[
  Log10@Norm[MatrixPower[A, n].P, Infinity], {n, 10, 1500},
  PlotRange -> {{10, Automatic}, All}
]

You are looking for the point where the above plot first crosses the axis, going from being >1 to being <1:

DiscretePlot[
  Log10@Norm[MatrixPower[A, n].P, Infinity], {n, 700, 870},
  PlotRange -> {{700, Automatic}, All}
]

It seems to me that you want to minimize the distance between your norm and $1$:

Clear[arg]
arg[n_?NumericQ] := (Norm[MatrixPower[A, n].P, Infinity] - 1)^2

DiscretePlot[Log10@arg[n], {n, 700, 870}, PlotRange -> {{700, Automatic}, All}]

NMinimize tries its best here, but unfortunately the minimum is not unique, so you do not necessarily get the lowest value of $n$ for which your condition is true:

NMinimize[
  {arg[n], n ∈ Integers, 750 < n < 810}, n,
  MaxIterations -> 450,
  Method -> "SimulatedAnnealing"
]

{0.000850499, {n -> 790}}

It seems that a manual search may indeed be a very good option here...

Answer 3

This is not as fast as MikeY's solution, but here it is.

A = SparseArray[{{1, 30} -> 1/2, {i_, j_} /; i - j == 1 -> 1}, {80, 80}] // N;
    P = Table[8 10^9 2^(-80) Binomial[80, k], {k, 0, 79}] // N;
    Reap[Do[If[Norm[MatrixPower[A, n].P, ∞] <= 1, Break[], 
 Sow[n + 1]], {n, 10000}]][[2, 1, -1]] // AbsoluteTiming

{0.146163, 801}

Or

n = 1;
While[True, If[Norm[MatrixPower[A, n].P, ∞] <= 1, Break[]]; n++]
n

801