How to make the computation of a minimum efficient?
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.
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.
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...
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