6
$\begingroup$

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 \[Element] 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 \[Element] Integers]

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

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – MarcoB Jan 13 at 22:07
  • $\begingroup$ @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. $\endgroup$ – Masacroso Jan 13 at 22:30
  • $\begingroup$ 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. $\endgroup$ – MarcoB Jan 13 at 22:41
10
$\begingroup$

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}]

enter image description here

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.

$\endgroup$
5
$\begingroup$

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}
]

log plot for integer values

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}
]

zoom in of above

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}]

Log arg

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...

$\endgroup$
  • $\begingroup$ 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 $\endgroup$ – Masacroso Jan 13 at 23:19
  • $\begingroup$ @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. $\endgroup$ – MarcoB Jan 13 at 23:23
3
$\begingroup$

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.