We can categorically state that we have not released man-eating badgers into the area
stečkin
S. B. Stečkin's lemma concerns behaviour of a trigonometric polynomial near its maximum. To be precise, if the polynomial $$ f(t) = \sum_{\abs{n}\leq N}\alpha_{n}e^{int} \quad (\alpha_{n}\in\complex) $$ is real-valued and $t_{0}\in [0,2\pi]$ satisfies $$ f(t_{0}) = \maxmod{f} = \max\left\{\abs{f(t)}\;:\; t\in [0, 2\pi]\right\}, $$ then $$ f(t_{0}+s) \geq \maxmod{f}\cos Ns \quad (\abs{s}\leq \pi/N). $$ In other words, the polynomial decays no faster that $\cos Ns$ close to a maximum without any reference to derivatives. This lemma and its generalisations can be used to construct “subdivide and reject” algorithms for global optimisation of non-convex objectives.
maxmod
The maxmod algorithm calculates the maximum modulus of a complex polynomial on the unit disc using a method of repeated subdivision of $[0,2\pi]$ and rejection of non-maximizing subintervals. A detailed description of the algorithm and its performance can be found here.
- maxmod.m for Octave/Matlab
- maxmod.py for Python
- a C99 implementation and its manual
- a fully templated C++11 implementation
The algorithm was used by Gregory Janks in a numerical demonstration that the Mandelbrot set is locally connected at the Feigenbaum point.
maxmodnd
Gabriel De La Chevrotière (McGill University) has generalised Stečkin's lemma to the multivariate case and so derives a method for finding the maximum modulus of a multivariate polynomial on the polydisc. The paper, published in SIURO, can be downloaded here.
- maxmodnd.m (requires polyvaln.m) for Octave/Matlab
- a version of the algorithm is included in the mvpoly package
maxmodnb
I have some experimental code for calculating the maximum modulus of a multivariate complex polynomial on the unit ball. Interested researchers are invited to contact me for a preliminary version of the code.
opnorm
S. W. Drury (McGill University) has proved a Stečkin-like estimate for operators from finite dimensional real vector spaces with the $p$-norm to Banach spaces with easily computed norms. This enables him to implement the first general algorithm for calculating the $\ell_p \rightarrow \ell_q$ operator-norms of matrices (albeit those with a small number of columns), and so allows him to find a counterexample to a long-standing conjecture of Matsaev.
- S. W. Drury's paper in Linear Algebra and its Applications (subscription required) and his webpage on the algorithm, including implementations for Visual C++ and Maple.
- The opnorm package, a multi-threaded C99 implementation (manual) targeted at Unix-like systems with bindings for Matlab, Octave and Python.
As can be seen from the table, the algorithm appears to require time which is exponential in the order of the matrix (in fact in the number of columns), so that only matrices with fewer that 10 columns are feasible. But this is to be expected, since Julien M. Hendrickx and Alex Olshevsky have shown (SIAM J. Matrix Anal. Appl. 01/2010; 31:2802–2812, DOI, arXiv) that the problem is NP-hard.
threads | ||||||||
---|---|---|---|---|---|---|---|---|
order | 1 | 2 | 4 | 8 | ||||
2 | 0.163 | 0.113 | 0.132 | 0.143 | ||||
3 | 0.949 | 0.530 | 0.349 | 0.383 | ||||
4 | 6.490 | 3.689 | 2.082 | 1.911 | ||||
5 | 49.962 | 27.154 | 15.061 | 12.445 | ||||
6 | 412.647 | 220.077 | 120.219 | 95.550 | ||||
7 | 3270.566 | 1746.786 | 946.903 | 740.649 |