# How to make the computation of a minimum efficient?

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.

• I don't fully understand your syntax or problem. Is it correct that you are looking for the minimum value of $n$ for which that Norm expression becomes less than 1? Or, similarly, for the maximum value of $n$ for which it is above 1? – MarcoB Jan 13 '20 at 22:07
• @MarcoB you are right, I wrote the code wrongly. My intention was to search for the minimal $n\in \Bbb N$ such that $\|A^nP\|_\infty <1$. Can you help me? I dont see how I can fix the code. – Masacroso Jan 13 '20 at 22:30
• Have you tried minimizing the Abs[Norm[ ... ] - 1], something like NMinimize[{Abs[Norm[ ... ] - 1], Element[n, Integers]}, n]? In other words, you would be trying to minimize the "distance of the norm from 1". I cannot test it myself at the moment unfortunately. – MarcoB Jan 13 '20 at 22:41

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[], P[], 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[]/2^(k - 1), P[]/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...

• thank you, but I want a way to calculate the minimal value of $n\in \Bbb N$, not of the norm, such that $\|A^n P\|_\infty<1$. This happens at $n=801$, that I can see giving values on the function g[n_]:=Max[MatrixPower[A,n].P]. This minimum is necessarily unique because the natural numbers are well-ordered – Masacroso Jan 13 '20 at 23:19
• @Masacroso I am not calculating a minimal value of the norm; I am calculating a minimal value of the distance of the norm from 1 though. Perhaps you might want to clarify what you seek in the OP, maybe I am misunderstanding. – MarcoB Jan 13 '20 at 23:23

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

• the advantage of this approach is that the function that works as a condition doesn't need to be monotone. Thank you. – Masacroso Jan 14 '20 at 15:07