main.tex 97.7 KB
Newer Older
Peter FRANEK's avatar
Peter FRANEK committed
%%%%%%%%%%%%%%%%%%%%% file template.tex %%%%%%%%%%%%%%%%%%%%%%%%%
Peter FRANEK's avatar
Peter FRANEK committed
2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134
% This is a general template file for the LaTeX package SVJour3
% for Springer journals.          Springer Heidelberg 2010/09/16
% Copy it to a new file with a new name and use it as the basis
% for your article. Delete % signs as needed.
% This template includes a few options for different layouts and
% content for various journals. Please consult a previous issue of
% your journal as needed.
% First comes an example EPS file -- just ignore it and
% proceed on the \documentclass line
% your LaTeX will extract the file if required
%!PS-Adobe-3.0 EPSF-3.0
%%BoundingBox: 19 19 221 221
%%CreationDate: Mon Sep 29 1997
%%Creator: programmed by hand (JK)
  20 20 moveto
  20 220 lineto
  220 220 lineto
  220 20 lineto
2 setlinewidth
  .4 setgray fill
%\documentclass{svjour3}                     % onecolumn (standard format)
%\documentclass[smallcondensed]{svjour3}     % onecolumn (ditto)
\documentclass[smallextended]{svjour3}       % onecolumn (second format)
%\documentclass[twocolumn]{svjour3}          % twocolumn
\smartqed  % flush right qed marks, e.g. at end of proof
% \usepackage{mathptmx}      % use Times fonts if available on your TeX system
% insert here the call for the packages your document requires
% etc.
\usepackage{courier} % this is for the algorithm listing...

% please place your own definitions here and don't use \def but
% \newcommand{}{}

%%%%%%%%%%%%%%%%%%% MACROS: Operators
\DeclareMathOperator{\im}{im} % IMAGE OF A MAP
\DeclareMathOperator{\id}{id}   % IDENTITY MAP
\newcommand\supp{{\mathop{\rm supp}\nolimits}}
\newcommand\sgn{{\mathop{\rm sgn}\nolimits}}
\newcommand{\argmax}{\mathop{\rm argmax}}
\newcommand{\conv}{\mathop{\rm conv}}
\newcommand{\aff}{\mathop{\rm aff}}


\def\d{\mathrm d}
\newcommand\makevec[1]{{\boldsymbol #1}}
\def \xx {\makevec{x}}
\def \yy {\makevec{y}}
\def \cc {\makevec{c}}
\def \ww {\makevec{w}}
\def \aa {\makevec{a}}
\def \mm {\makevec{m}}
\def \gg {\makevec{g}}
\def \nn {\makevec{n}}
\def \nula {\makevec{0}}
\def \ll {\boldsymbol{\ell}}

\newcommand{\heading}[1]{\vspace{1ex}\par\noindent{\bf\boldmath #1}}

Peter FRANEK's avatar
Peter FRANEK committed
135 136
Peter FRANEK's avatar
Peter FRANEK committed
137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173
\vspace{1mm}#1 \vspace{1mm}

\def\PSPers{\textbf{Primary--Secondary Persistence}}
\def\PersPrim{{\textbf{ObstructionPersistence}$_{\boldsymbol 1}$}}
\def\PersSec{{\textbf{ObstructionPersistence}$_{\boldsymbol 2}$}}
\def\Obst#1{{\textbf{Obstruction}$_{\boldsymbol #1}$}}


\newcommand{\cmt}[1]{{\color{Orange} #1}}
{\color{Orange}\textsf{*** (MAREK: ) #1\\}}}

{\color{blue}\textsf{*** (PETER: ) #1\\}}}
\newcommand{\HW}[1]         {\marginpar{\textcolor{black}{hubert}} {{\textcolor{blue}{#1}}}}


% Insert the name of "your journal" with
% \journalname{myjournal}

\title{Solving Equations and Optimization Problems with Uncertainty\thanks{The research 
Peter FRANEK's avatar
Peter FRANEK committed
175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198
leading to these results has received funding from Austrian Science Fund (FWF): M 1980 and
from the Czech Science Foundation (GACR) grant number 15-14484S with institutional support RVO:67985807.}
%\subtitle{Do you have a subtitle?\\ If so, write it here}

%\titlerunning{Short form of title}        % if too long for running head

\author{Peter Franek \and
        Marek Kr\v c\'al \and \hbox{Hubert Wagner}

%\authorrunning{Short form of author list} % if too long for running head

\institute{IST Austria

\date{Received: date / Accepted: date}
% The correct dates will be entered by the editor


We study the problem of detecting zeros of continuous functions that are known only up to an error bound, extending
Peter FRANEK's avatar
Peter FRANEK committed
the theoretical work of~\cite{nondec} with explicit algorithms and experiments with an implementation.\footnote{\url{}}
Peter FRANEK's avatar
Peter FRANEK committed
200 201
%More formally, the \emph{robustness of zero} of a continuous map $f\:X\to \R^n$ is the maximal $r>0$ such that  each $g\:X\to\R^n$ 
%with $\|f-g\|_\infty\le r$ has a zero. We develop and implement an algorithm approximating the robustness of zero. 
Peter FRANEK's avatar
Peter FRANEK committed
202 203 204 205
%The main source of motivation for such an algorithm comes from the field of interval arithmetic.
Further, we show how to use the algorithm for approximating worst-case optima in optimization problems in which 
the feasible domain is defined by the zero set of a function $f: X\to\R^n$ which is only known approximately.

Peter FRANEK's avatar
Peter FRANEK committed
206 207 208 209 210 211 212
The algorithm first identifies a~subdomain $A$ where the function $f$ is provably non-zero, 
a~simplicial approximation $f': A\to S^{n-1}$ of $f/|f|$, and then verifies non-extendability of $f'$ to $X$ to certify a~zero.
Deciding extendability is based on computing the cohomological
\emph{obstructions} and their persistence. We describe an explicit algorithm for the \emph{primary and secondary 
obstruction}, two stages of a sequence of algorithms with increasing complexity. Using elements and techniques of persistent homology,
we quantify the persitence of these obstructions and hence of the robustness of zero.

Peter FRANEK's avatar
Peter FRANEK committed
213 214 215 216 217 218 219 220 221 222 223 224 225 226 227
We provide experimental evidence that for random Gaussian fields,
the \emph{primary obstruction}---a much less computationally demanding test than the secondary obstruction---is typically 
sufficient for approximating robustness of zero.

\keywords{Computational homotopy theory \and Uncertainty \and Obstructions }
% \PACS{PACS code1 \and PACS code2 \and more}
% \subclass{MSC code1 \and MSC code2 \and more}
Detecting zeros of $\R^n$-valued functions is  equivalent to solving systems of real equations, a fundamental problem
of mathematics and theoretical computer science. Our research is  motivated by practical applications, in which the data is often known only approximately. 
Peter FRANEK's avatar
Peter FRANEK committed
We address the case where the input data is limited to the  \emph{approximate} values of a~continuous function $f$ with
Peter FRANEK's avatar
Peter FRANEK committed
values in $\R^n$, sampled over a finite point set.
Peter FRANEK's avatar
Peter FRANEK committed
230 231
This uncertainty  is handled in a deterministic way: we aim at verifying that \emph{each} continuous
function compatible with our partial knowledge of $f$, has a~zero. 
Peter FRANEK's avatar
Peter FRANEK committed
232 233

Functions that are known only approximately appear in various contexts and are handled in different ways. 
Peter FRANEK's avatar
Peter FRANEK committed
For example, rounding errors in floating-point computations are systematically treated by methods of \emph{interval}
Peter FRANEK's avatar
Peter FRANEK committed
235 236
Peter FRANEK's avatar
Peter FRANEK committed
237 238
\caption{For scalar valued function, existence of a~zero can be verified via the intermediate value theorem.}
Peter FRANEK's avatar
Peter FRANEK committed
Peter FRANEK's avatar
Peter FRANEK committed
\emph{arithmetic} and detection of zeros 
Peter FRANEK's avatar
Peter FRANEK committed
241 242 243 244
resistant to bounded errors is a~frequent problem in this field \cite{Frommer:05,Dian:03,Alefeld:04,Jaulin1,Jaulin4}. 
Other instances of uncertain functions come from measurements of physical quantities, such as in medical 
or robotics~\cite{Merlet:2009,Jaulin3}.
Peter FRANEK's avatar
Peter FRANEK committed
We suppose that potential applications include robust detection of level sets $f^{-1}(a)$ in medical image processing, 
Peter FRANEK's avatar
Peter FRANEK committed
246 247 248 249 250
analysing robot trajectory based on data obtained from sensors~\cite{Jaulin3,Jaulin4},
or computing the \emph{inner approximation} of reachable regions of a robotic arm~\cite{Jaulin5}.
The algorithm could also be exploited for analysis of functions obtained by regression (say, in machine learning),
where the function is chosen to fit some given set of sampled values.

Peter FRANEK's avatar
Peter FRANEK committed
251 252 253
To verify that a function $f$, of which we only have a~limited knowledge, has a zero, is equivalent to showing that \emph{each} 
potential candidate $g$ for $f$ has a zero.  If we only have access to sampled values of $f$ and a~Lipschitz constant,
then the set of all such admissible functions $g$ is huge and can not be finitely parametrized. 
Peter FRANEK's avatar
Peter FRANEK committed
254 255
However, methods of computational homotopy theory can be applied: 
the closely related problem of verifying that each continuous 
Peter FRANEK's avatar
Peter FRANEK committed
256 257 258 259 260 261
$r$-perturbation of a~given function has a~zero, can be reduced to the topological extension problem for maps into a sphere~\cite{nondec}. 
The latter problem can be addressed via means of obstruction theory, using an algorithmic construction of Postnikov towers.
Such construction have never been implemented and is in its full generality probably out of reach, given the limitations of
computer power. Using a number of simplifications as well as some methods of persistent homology, we present
a~partial solution to the above problem accompanied by an implementation, complexity analysis and several computational experiments.

Peter FRANEK's avatar
Peter FRANEK committed

Peter FRANEK's avatar
Peter FRANEK committed
263 264 265 266 267 268 269 270 271 272 273
%--- ? ---
%The generality of our result, combined with its algorithmic efficiency, open doors for new applications. 
%One of them is for optimization problems: given an objective function $u$ and a partial knowledge of $f: X\to\R^n$, 
%we can compute a lower bound for $\max \, u(x)$ subject to $f(x)=0$. 
%The best such lower bound is called~\emph{robust counterpart}: we can approximate it in cases when $f$ is given via function values in a finite set 
%and a~Lipschitz constant.
%%More broadly, we aim to present computational homotopy theory as a useful tool and a fascinating area of 
%%algorithmic research to a~wide computer science audience. 
Peter FRANEK's avatar
Peter FRANEK committed
274 275 276 277 278 279 280 281 282 283 284 285 286 287 288
\heading{Statement of the results.} 
We present an algorithm for detecting zeros of vector valued functions $f:X\to\R^n$ on a compact space $X$ and for approximating the 
\emph{robustness} of zero, that is, a~maximal real number $r>0$ such that every continuous $g: X\to\R^n$ satisfying $\|g-f\|\leq r$
has a zero. By $\|f\|$ we denote the max norm $\max_{x\in X} |f(x)|$ where $|\cdot|$ is a fixed $\ell_p$ norm in $\R^n$. 
Nontrivial cases happen if $\dim X\geq n$, as otherwise arbitrarily small perturbations of $f$ avoid zero. 
For computer representation we assume that the space $X$ is a simplicial complex.
The map $f\:X\to\R^n$ is specified by its values on the vertices which are assumed to be rational,
and by a rational value $\alpha>0$ such that $|f(x)-f(y)|\le \alpha$ for arbitrary points $x$ and $y$ of any simplex of $X$. 
%Such a function will be referred to as \emph{simplexwise $\alpha$-Lipschitz.} 
We emphasize that the precise knowledge of $f$ is not needed.
The algorithm computes two numbers $r,R$ such that
\item Every continuous $g$, $\|g-f\|\leq r$, has a zero, and
\item Some continuous $g$, $\|g-f\|\leq R$, has no zero.
Peter FRANEK's avatar
Peter FRANEK committed
289 290
We will say that $f$ has an \emph{$r$-robust zero}.
A positive $r>0$ is then of course a~certificate of existence of zero of $f$.
Peter FRANEK's avatar
Peter FRANEK committed
291 292 293 294
In the dimension range $\dim X\leq n+1$ or $n<3$ the gap $R-r$ provably converges to zero, 
if the constant $\alpha$ (and hence our lack of knowledge of $f$) goes to zero. However, the algorithm outputs a~correct
and reasonable lower bound $r$ of robustness for inputs of arbitrary dimensions.

Peter FRANEK's avatar
Peter FRANEK committed
295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312
Our second result is based on computational experiments with random functions. A natural informal question is
\item {\it
How typical are functions for which higher obstructions are needed for detecting a zero?}
An example of a~function with nontrivial secondary obstruction is any maps $f: B^4\to\R^3$ such that $f|_{\partial B^4}$
is homotopic to the Hopf fibration $S^3\to S^2$ and hence cannot be extended to $B^4\to \R^3\setminus \{0\}$.
Such property can be verified, if we are given a~sample of function values and a Lipschitz 
constant. Moreover, the homotopy class of $f|_{\partial B^4}$ does not change if we slightly perturb these sampled function values.

Surprisingly, when performing experiments with random functions (mainly random Gaussian fields), we observed that higher obstructions
are typically \emph{not} needed. 
Whenever we detected a zero of a~randomly generated uncertain function, it was via means of primary obstruction only.
We performed experiments with various random functions from a~triangulated $4$-cube or from a~$4$-torus 
into $\R^3$ as well as from a $5$-torus into $\R^4$.
This observation---if confirmed by theory or by more experiments in different 
Peter FRANEK's avatar
Peter FRANEK committed
313 314 315 316 317 318 319 320 321 322 323 324 325 326 327
settings---could justify the usage of \emph{only the primary obstruction} in potential 
future engineering applications. 

\heading{State of the art.}
Algorithms for detecting zeros used in software packages are based on iterative methods which are often applicable 
if $f$ is given by formulas and is differentiable.
However, these algorithms usually give no guarantees of correctness: the satisfiability of $f(x)=0$ is undecidable for any
class of real functions $f$ that contain polynomials and the sine function~\cite{Wang:74}.

A number of methods has been proposed for testing the (non-)existence of zeros of continuous functions,
exploiting tools ranging from iterative methods in numerical analysis to topology.
The problem has been most studied in the case $\dim X=n$. If  $B^n$ is a unit ball in $\R^n$, then
verifying zeros of $f: B^n\to\R^n$ is equivalent to verifying a fixed point of $f+\mathrm{id}$: here the Brouwer fixed point theorem can be applied~\cite{Rump:2010}.
Other methods for zero verification were studied in the field of interval arithmetic, such as \emph{Miranda's test}~\cite{Alefeld:01},
\emph{Borsuk's test}~\cite{Frommer:05} and the \emph{degree test}~\cite{Franek-al}.
Peter FRANEK's avatar
Peter FRANEK committed
328 329
All of these tests have topological flavour and are stable with respect to perturbations of the input function. 
Authors of~\cite{Franek_Ratschan:2015} show that the \emph{degree test} can detect a zero of $f: B^n\to\R^n$ whenever
Peter FRANEK's avatar
Peter FRANEK committed
330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356
the zero of $f$ is \emph{robust} (that is, each $g$ close enough to $f$ has a zero). 
The above mentioned primary obstruction directly reduces to the degree test if $\dim X=n$.
While the topological degree computation has been explicitly described in the literature and also implemented~\cite{degpaper}, 
the problem is far more complicated if the domain $X$ has larger dimension than $n$.
%(then the zero set is the solution of an~underdetermined system of equations). 
In~\cite{nondec} we showed that the existence of a~robust zero of piecewise linear functions $X\to\R^n$ is undecidable 
if $n\geq 3$ is a fixed odd integer and the input $X$ is an~arbitrary $(2n-2)$-dimensional simplicial complex.

Zero sets of functions with inherent uncertainty were studied via means of computational topology 
in the context of~\emph{well groups}~\cite{ComputeWell}.
In their general settings, well groups associated to $f: X\to Y$ and a subspace $Y'\subseteq Y$ 
describe properties of the preimage $f^{-1}(Y')$ which persist if we perturbate the input function $f$. 
In the important case of $Y=\R^n$ and $Y'=\{0\}$, well groups describe zero sets of functions: namely, the zeroth well
groups measures robustness of \emph{existence} of zero and higher well groups reflect further topological properties of the zero sets.
In~\cite{well-socg} we have shown that the primary obstruction can be used to compute certain subgroup of the well group
that in many cases coincides with the full well group (see~\cite[Thm. 1.4]{well-socg}).
A general algorithm for well group computation cannot be expected because the above mentioned undecidability result~\cite{nondec}
directly transfers  to well groups when $\dim X\ge 2n-2$ as well. 
Our implementation can be thought of as~an approximation of the zeroth well group that extends the work of~\cite{chazal} 
where the special case $\dim X=n$ is solved. 

An application of our algorithm is in worst-case analysis of optimization problems where the feasible domain is defined by equations.
Worst-case approach in robust optimization has been widely studied, see~\cite{Ben-Tal2002,Ben:2009,bertsimas2011theory,Beyer20073190}. 
Usually, the uncertainty applies to a finite number of parameters which are assumed to be taken from a~known domain. 
In our approach, we rather work with the space of \emph{all continuous functions} that are compatible
with our partial knowledge of $f$. %: the tradeoff is that we don't aim at instances of $f: X\to\R^n$ where $\dim X$ is high.

Peter FRANEK's avatar
Peter FRANEK committed
\heading{Outline and organization of the paper.}
Peter FRANEK's avatar
Peter FRANEK committed
358 359 360 361 362
The main tool is to find a domain $A\subseteq X$ where $f$ is provable nonzero and where our knowledge of $f$ is sufficient
to determine the homotopy class of $f|_A$ as a~map to $\R^n\setminus \{0\}\simeq S^{n-1}$.
Then non-extendability to $X\to S^{n-1}$ is a certificate of zero.

Peter FRANEK's avatar
Peter FRANEK committed
In the algorithm, we first create a~finite filtration $\{A_r\supseteq A_s\}_{r\leq s}$ of subcomplexes that ``approximate'' the topological spaces 
Peter FRANEK's avatar
Peter FRANEK committed
364 365 366 367
$A(r):=\{x\in X:\,\,|f(x)|\geq r\}$. 
We compute a~simplicial approximation $f': A_r\to \Sigma$ of $f$ where $\Sigma$ is a given triangulation
of the $(n-1)$-sphere.
Then we ask for the smallest $r$ such that the restriction of $f'$ to $A_r$ can be extended to
all of $X$, and show that the robustness of zero of the original function is $\alpha$-far from $r$.
Peter FRANEK's avatar
Peter FRANEK committed
369 370 371 372 373 374 375 376

Such extendability is decidable if $\dim X\leq 2n-3$~\cite{ext-hard}, but the only procedere for this we are aware of
is based on the~algorithm for computing stages of Postnikov towers from~\cite{polypost} 
that depends on several other papers~\cite{pKZ1,post,vokrinek:oddspheres} and is unlikely to be fully implemented in near future. 
Instead of that, we implemented a~persistent version of both the \emph{primary} and \emph{secondary obstruction}, 
which test extendability to the $n$- and $(n+1)$-skeleton of $X$.\footnote{The only exception is the case $n=3$, $\dim X>3$ where
the triviality of secondary obstruction is undecidable in general. However, if $X$ is assumed to be a triangulation of the cube $[0,1]^4$,
then our algorithm works with no essential changes. For many other fixed $4$-dimensional spaces $X$ the problem is decidable too.} 
Peter FRANEK's avatar
Peter FRANEK committed
First, we compute the minimal $r_1$ for which the cohomological obstructions to extending 
Peter FRANEK's avatar
Peter FRANEK committed
${f'}|_{A_{r_1}}$ to $A_{r_1}\cup X^{(n)}$ (\emph{primary obstruction}) vanishes.
Peter FRANEK's avatar
Peter FRANEK committed
Similarly, we compute a~minimal $r_2\geq r_1$ for which $f'|_{A_{r_2}}$ is extendable to $A_{r_2}\cup X^{(n+1)}$ 
Peter FRANEK's avatar
Peter FRANEK committed
380 381 382 383
(vanishing of the \emph{secondary obstruction}): this requires to parametrize all extensions to the $n$-skeleton.
%Via computational experiments with random Gaussian fields in low dimensions, we observe that 
%the \emph{primary obstruction}, measuring non-extendability to the $n$-skeleton, is typically sufficient to detect a~robust zero of $f$.\footnote{This
%observation was somehow surprising, given that the set of functions for which $r_2>r_1$, is open in $C(X,\R^n)$.}
Peter FRANEK's avatar
Peter FRANEK committed
384 385 386 387

In Section~\ref{s:discretize} we show how to approximate the spaces $A(r)$ and sphere-valued maps $f/|f|$ via simplicial functions. 
A~high-level description of our algorithm is in Sections~\ref{sec:algorithm-oracle} and \ref{s:obstructions},
with a partial~lower-level description in Appendix~\ref{s:secondary-persistence} and \ref{s:implementation}.
Peter FRANEK's avatar
Peter FRANEK committed
388 389
In Section~\ref{s:rob-opt} we show how to use the method for approximating
worst-case optima in optimization problems where the feasible domain is defined by equations.
Peter FRANEK's avatar
Peter FRANEK committed
390 391
In Section~\ref{s:experiment} we present some computational experiments with random Gaussian fields.
More details about testing and performance are delegated to Appendix~\ref{s:formulas}. 
Peter FRANEK's avatar
Peter FRANEK committed
The last section contains a~theoretical worst-case complexity bounds.
Peter FRANEK's avatar
Peter FRANEK committed
393 394 395

\section{Discretizing the function $f$.}
Peter FRANEK's avatar
Peter FRANEK committed
Peter FRANEK's avatar
Peter FRANEK committed
397 398 399 400 401
In this section we show how to convert the ``unknown'' continuous function $f$ to its discrete simplicial approximation.
A \indef{continuous filtration} of spaces is a family $(A_r)_{r\in\R}$ such that $A_r\supseteq A_s$ whenever $r\le s$.
Peter FRANEK's avatar
Peter FRANEK committed
A continuous filtration $(A_r)_{r\in\R}$ is called \indef{step-like} whenever there exists a sequence of numbers $-\infty=:r_{-1}<r_0\le r_1\le r_2\le\ldots\le r_k$ such that for any $r,s\in (r_i,r_{i+1}]$, $A_r=A_s$ holds for all $i$. 
Peter FRANEK's avatar
Peter FRANEK committed
403 404 405 406 407 408 409 410
Continuous filtrations $(A_r)_r$ and $(B_r)_r$ are called \indef{$\alpha$-interleaved} whenever $B_{r+\alpha}\subseteq A_r$ and $A_{r+\alpha}\subseteq B_r$ for each $r\in\R$.

Let $f\:X\to \R^n$ be a continuous map on a simplicial complex $X$ and let $|\cdot|$ be a norm on $\R^n$. 
Peter FRANEK's avatar
Peter FRANEK committed
\item By $A_r$ we denote the \emph{subcomplex} of $X$ spanned by the vertices $v$
Peter FRANEK's avatar
Peter FRANEK committed
of $X$ with $|f(v)|\ge r$. 
Peter FRANEK's avatar
Peter FRANEK committed
\item By $A(r)$ we denote the \emph{subspace} of $X$ defined by $A(r)=\{x\in X\: |f(x)|\ge r\}$. 
Peter FRANEK's avatar
Peter FRANEK committed
414 415 416 417 418
\item We say that $f$ is \indef{simplexwise $\alpha$-Lipschitz} whenever $|f(x)-f(y)|\le \alpha$ for each pair of points $x,y\in\Delta$ of any simplex $\Delta\in X$.

Spaces $A_r$ form a step-like filtration where the steps occur for each $r$ equal to $|f(v)|$ for some vertex $v$ of $X$.
Peter FRANEK's avatar
Peter FRANEK committed
Let $\Sigma^{n-1}$ be the simplicial model of the $(n-1)$-sphere obtained from the boundary of a cross-polytope.
Peter FRANEK's avatar
Peter FRANEK committed
420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484
%Simplices of $\Sigma^{n-1}$ are all subsets of $\{e_1,-e_1\ldots,e_n,-e_n\}$ that do not contain a~pair of antipodals $\pm e_j$.
A natural sphere-valued approximation of $f$ is then given by the map $f'$ as follows.
Let $f: X\to\R^n$ and $V$ be a subset of the vertices of $X$. We define the \indef{vertex approximation} 
$f': V\to \{e_1, -e_1,\ldots, e_n, -e_n\}$ to by the map that to a vertex $v$ assigns $s_j e_j$, where
$j$ is the index of the component of $f(v)$ with \emph{largest} absolute value and $s_j$ is the sign of $f_j(v)$.\footnote{For example,
if $f(v)=[2,-3]$, then we choose $f'(v)=-e_2$. If there are more components of $f(v)$ with the same absolute value, we choose one by arbitrary
Let $f\:X\to \R^n$ be a simplexwise $\alpha$-Lipschitz map for some constant $\alpha>0$, $p\in [1,\infty]$ and $|\cdot |$ be the $\ell_p$-norm. 
Then the following holds: 
\item The continuous filtrations $(A_r)_{r\in \R}$ and \(A(r)_{r\in\R}\) are $\alpha$-interleaved. 
\item If $r> \alpha n^{1/p}/2,$ the vertex approximation $f': V(A_r)\to V(\Sigma^{n-1})$
defines a simplicial map $f'\:A_r\to\Sigma^{n-1}$ (that is, it maps simplices to simplices). 
\item If $r>\alpha n^{1/p}$, then $f'\:A_r\to \Sigma^{n-1}\subseteq \R^n\setminus\{0\}$ 
is homotopic to $f|_{A_r}\:A_r\to \R^n\setminus\{0\}$.

% Motivated by one of SoCG reviewers.. :
The simplicial map $f'\:A_r\to\Sigma^{n-1}$ as above will be called the \indef{simplicial approximation of $f|_{A_r}$.}
We note that the choice of the sphere model and the discretization $f'$ of $f$ is independent of the rest of the algorithm given in the following chapters
and is not the only possible choice: however, we couldn't find one with approximative properties better than in Lemma~\ref{t:approx}.
\begin{proof}[Proof of part 1]
If a~point $x\in |X|$ is not in $A_r$, we know that  $|f(v)|<r$ for some vertex $v$ of a~simplex $\Delta$ supporting $x$. 
The simplexwise $\alpha$-Lipschitz property then implies that $|f(x)|<r+\alpha$, hence $x\notin A(r+\alpha)$. This proves
$A(r+\alpha)\subseteq A_r$.

The other inclusion holds because once a point $x$ of a simplex $\Delta\in X$ is in a~given $A_r$, 
then for arbitrary vertex $v$ of $\Delta$ holds $|f(v)|\ge r$ and  $|f(x)-f(v)|\le \alpha$. 
Thus $|f(x)|\ge r-\alpha$, hence $x\in A(r-\alpha)$. 
\begin{proof}[Proof of part 2]
We want to prove that no adjacent vertices $u$ and $v$ are mapped by $f'$ to $e_i,-e_i$ for some $i$. 
Without loss of generality we can assume that $f'(v)=e_1$. Assuming $|f(v)|\geq r$ and the definition of $f'$, 
the first component (with largest absolute value) must satisfy $f_1(v)\ge rn^{-1/p}$.
Similarly if $f'(u)=-e_1,$ then $f(u)_1\le-rn^{-1/p}$, but this would imply that $\alpha\geq |f(u)-f(v)|\ge2rn^{-1/p}$, 
contradicting the assumption $r>\alpha n^{1/p}/2$.
\begin{proof}[Proof of part 3]
We show that the simplexwise straight-line homotopy between $f$ and $f'$ has values in $\R^n\setminus \{0\}$. 
Let $\Delta\in X$, $v$ be a vertex of $\Delta$, $x\in\Delta$, and assume that $r>\alpha n^{1/p}$. Again, assume WLOG that $f'(v)=e_1$.
Then $f_1(x)\geq f_1(v) - \alpha \geq r n^{-1/p}-\alpha >0$ and the straightline homotopy between maps 
$f,f'\:A_r\to\R^n\setminus\{0\}$ has the first coordinate positive, hence it avoids zero. 

\marek{Nasledujici jsem psal pro pripad simplexwise linear. Pro ted si myslim, ze to plati i v nasem soucasnem obecnejsim pripade, ale jeste double check.}We remark that the weaker assumption $r>\alpha
n^{1/p}/2$ cannot be further weakened: for instance
when $p=2$ and when a pair of adjacent vertices are mapped
to $r/\sqrt n(1+\epsilon,\ldots,1)$ and $r/\sqrt n(-1-\epsilon,1,\ldots,1)$, then their simplicial assignment is $e_1$ and $-e_1$ for any $\epsilon>0$. The
stronger assumption can be weakened. In the case of the $\ell_2$
norm, we can show that the assumption $$r>{\sqrt{n}\over 2\sqrt{n-1}} \alpha \sqrt{n}$$ is sufficient. But the best example we know where the simplicial map $f'$ is
correctly defined but not homotopic to $f$ on $A_r$ achieves only $$r={1\over 2\sqrt{2}}
\alpha\sqrt{n}.$$ We conjecture that this is the optimal bound.
Proving this optimal bound would be a practically relevant achievement since
the simpliciality can be checked algorithmically relatively fast and in our instances holds  for lower values of $r$ very often. On the other hand, checking the homotopy correctness would be computationally demanding so we have to rely on this theoretical bound if we want to guarantee the correctness of our algorithm.

\section{The algorithm using an oracle for persistence of obstructions.}
Peter FRANEK's avatar
Peter FRANEK committed
In this section we describe a high-level description of our algorithm for approximating robustness of zero. It uses an oracle for computing the 
Peter FRANEK's avatar
Peter FRANEK committed
486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672
quantities introduced the next definition.
Let $X\supseteq A_{r_0}\supseteq A_{r_1}\supseteq \ldots$ be a filtration of simplicial complexes, 
$r_i\leq r_j$ for $i\leq j$ and $f': A_{r_0}\to S^{n-1}$. 
Then the \indef{persistence of primary obstruction} is the smallest $r_j$ such that the restriction of $f'$ to $A_{r_j}$ 
\emph{can} be extended to a~(not necessarily simplicial) map $A_{r_j}\cup X^{(n)}\to S^{n-1}$ where $X^{(n)}$ is the $n$-skeleton of $X$. The \indef{persistence of secondary obstruction} is the smallest $r_k$
such that the restriction of $f'$ to $A_{r_k}$ can be extended to $A_{r_k}\cup X^{(n+1)}\to S^{n-1}$.
In what follows, assume that an oracle is given that, for a filtration of simplicial complexes and a simplicial map $f'$, 
computes the persistence of secondary obstrucion. 
We assume that we are given a continuous map $f\:X\to\R^n$ by its values on the vertices of $X$,
its simplexwise Lipschitz constant $\alpha$ and a norm $\ell_p$ on $\R^n$ for $p\in [1,\infty]$.
%\footnote{If we run the algorithm 
%with the parameter $p=1$, it the outputs are correct for arbitrary norm satisfying $|e_i|=1$ for all $i$. 
%However, for larger values of $p$ it gives better approximation guarantees.} 
The outline of the algorithm follows.

%Compute the edgewise Lipschitz constant $\alpha = \max_{(v,u)\text{ %edge of } X} |f(u)-f(v)|$ of $f$.
\item \begin{enumerate}
Label the set $\{|f(v)|\:v\in V(X)\text{ such that } |f(v)|\geq \alpha n^{1/p} \}$   by $\{r_0,r_1,\ldots,r_h\}$ 
so that $r_i\leq r_j$ for $i\leq j$.
For any simplex $\Delta\in X$ compute its filtration value $r(\Delta)$ by 
r(\Delta):=\min_{v\text{ vertex of }\Delta} |f(v)|.
This yields a filtration $A_{r_0}\supseteq\ldots\supseteq A_{r_h}$ that together with the values $r_0,\ldots,r_h$ 
determines the step-like continuous filtration $(A_r)_{r\in\R}$ from Definition~\ref{d:filtration}.
For vertices $v$ of $X$ with $|f(v)|\ge r_0$ compute the vertex approximation $f'(v)$ via Definition~\ref{d:vertex-app}. 
\item Use the oracle to compute the persistence of secondary obstruction $r_k$ (Def.~\ref{d:persistence-obstruction}). 
\item {\bf If} $k>0$: output ``robustness of zero is at least  $r_{k}-\alpha.$''\\
      {\bf Else}: output ``No guarantee of zero.''
\item {\bf If} $\dim X\leq n+1$ or $n<3$:  output ``robustness of zero is at most $r_k+\alpha$'' \\
      {\bf Else}: output ``No guarantee of upper bound of robustness''

The constraints in {\bf  B.b.} could be replaced by $\dim X\leq n-1+k$, if we used an oracle for persistence of the first $k$ obstructions.
However, implementing such oracle is theoretically possible only if $k<n-1$. In fact, even the special case $\dim X=4$ and $n=3$
is beyond this bound and we cannot implement the oracle for secondary obstruction for this dimension pair
with no restrictions on $X$. In the important special case when $X$ is topologically a cube $[0,1]^4$ and $n=3$, the general algorithm works 
with no essential changes. More details about the implementation of the oracle for this dimension pair
are given in Appendix~\ref{s:secondary-persistence}, p.~\pageref{h:mn=43}.

\begin{theorem}\label{t:main} The above algorithm outputs correct statements. 

In~\cite[Lemma 3.3]{nondec} we showed that $f$ has an $r$-robust zero iff $f|_{A(r)}$ is not extendable to a nowhere zero function on $X$.

{\it Correctness of {\bf B.a.}}
Assume first that $k>0$. That means that $f'|_{A_r}$ is not extendable to a~function $X\to\R\setminus\{0\}$ 
for $r<r_k$. For any $r\in (r_0, r_k)$, $f'|_{A_r}$ is well defined and homotopic to $f|_{A_r}: A_r\to\R\setminus \{0\}$ 
by Lemma~\ref{t:approx}, 3. It follows that $f|_{A_r}$ is not extendable to a nowhere zero function $X\to\R^n\setminus\{0\}$ either.
Further, by Lemma~\ref{t:approx}, 1, $A_r\subseteq A(r-\alpha)$ and the non-extendability of $f|_{A_r}$ immediately implies
the non-extendability of $f|_{A(r-\alpha)}$. Thus we conclude that $f$ has at least $(r-\alpha)$-robust zero for arbitrary $r\in (r_0, r_k)$.
The set $\{r>0 \,\, | \,\, \,\text{each }\,\, g,\,\, \|g-f\|\leq r,\,\text{ has a zero}\}$ is closed (see~\cite[p. 9]{nondec} for an argument)
and hence $f$ has at least $(r_k-\alpha)$-robust zero.

{\it Correctness of {\bf B.b.}}
The restriction of $f'$ to $A_{r_k}$ is extendable to $A\cup X^{(n+1)}$, so for any $r\geq r_k$, 
$f'|_{A_{r}}$ is extendable to $A_r\cup X^{(n+1)}$ as well. If $\dim X\leq n+1$, then this is equivalent to
the extendability to all of $X$. 
The cases $n<3$ reflect low dimensional phenomena: we will show that then the extendability to $A\cup X^{(n)}$ implies the extendability 
to all of $X$.
If $n=1$, $f'$ has values in the $0$-sphere $S^0\in \{+,-\}$ and if it 
can be extended to the $1$-skeleton, we can assign a sign $+$ or $-$ to each connected component of $X$ and naturally extend to $X\to \{+,-\}$. 
If $n=2$, then $f'$ has values in a circle, $S^1$. Assume that it can be extended to $A\cup X^{(2)}\to S^1$.
Then any extension $A\cup X^{(j)}\to S^1$ of $f'$, $j\geq 2$, can be extended to $A\cup X^{(j+1)}$, because the restriction 
of $g$ to the boundary of any $(j+1)$-simplex $\Delta^{j+1}$ defines a map $\partial \Delta^{j+1}\to S^1$
from a $j$-sphere to the circle and such map is homotopic to a constant (see, e.g.~\cite[Chapter 4.1]{Hatcher}),
hence $\partial \Delta^{j+1}\to S^1$ can always be extended to all of $\Delta^{j+1}$.

Assume that $\dim X\leq n+1$ or $n<3$ and that $r>r_k$. 
Then $r>r_0$ and $f'|_{A_r}$ is well defined and homotopic to $f|_{A_r}$ by Lemma~\ref{t:approx}.
This implies the extendability of $f|_{A_r}$ to all of $X$ 
and the relation $A(r+\alpha)\subseteq A_r$ implies the extendability of $f|_{A(r+\alpha)}$. Thus the robustness is less than $r+\alpha$ for any
$r>r_k$, yielding an upper bound $r_k+\alpha$ on the robustness of zero.

\section{Persistence of obstructions}
In this section we describe the algorithm for computing the persistence of primary obstruction and roughly outline the algorithm for secondary 
obstruction (which is described in more detail in Appendix~\ref{s:secondary-persistence}).

\heading{Primary obstruction---extendability to $X^{(n)}$.} 
Here we review some facts from obstruction theory. 
A reference for the next proposition can be a~textbook such as \cite[III, 1.2]{prasolov}. 
\begin{proposition}[Primary obstruction]\label{p:prim}
Let $A\subseteq X$ be a pair of simplicial complexes, $f\:A\to \Sigma^{n-1}$ be simplicial, $z\in C^{n-1}(\Sigma^{n-1};\Z)$ be a cocycle
generating the cohomology, and $y:=f^ \sharp(z) \in Z^{n-1}(A;\Z)$ its pullback. 

Then $f: A\to \Sigma^{n-1}$ can be extended to a (not necessarily simplicial) map $A\cup X^{(n)}\to\Sigma^{n-1}$, iff 
$y\in Z^{n-1}(A;\Z)$ can be extended to a~cocycle $x\in Z^{n-1}(X;\Z)$ such that $x|_A=y$.

Thus extendability of $f$ is reduced to extendability of an $A$-cocycle $y$ to a global cocycle $x$ defined on all of $X$.
We will use the notation $\Omega(A):=\{x\in Z^{n-1}(X;\Z)\: x|_A = y\}$ of all cocycle extensions and want to test its non-emptiness.
$\Omega(A)$ corresponds to solutions of a~linear equation over integers. 
To see that, let $\bar y\in C^{n-1}(X;\Z)$ be an arbitrary cochain (not necessarily a cocycle) such that $\bar y|_A=y$. We have that
\Omega(A)=\{\bar y-c\: c\in C^{n-1}(X,A;\Z)\text{ such that }\delta c = \delta \bar y\}.\end{equation} 
Subtracting $c\in C^{n-1}(X,A;\Z)$ does not change the values on $A$-simplices, so any such $\bar y-c$ is still an extension of $y$.
The non-emptiness of $\Omega(A)$ is thus equivalent to solvability of the linear equation  $\delta c=\delta \bar y$ with the 
unknown $c\in C^{n-1}(X,A;\Z)$. 

A natural set of generators of $C^{n-1}(X,A;\Z)$ is the set of all $(n-1)$-simplices in $X$ that are not in $A$
with the identification between a simplex $\Delta$ and its \emph{characteristic cochain} that assigns $1$ to $\Delta$ 
and $0$ to all other simplices.  Converting $\delta c=\delta\bar y$ into an explicit matrix system of linear equations 
then amounts to enumerating the $(n-1)$- and $n$-simplices in $X\setminus A$, computing the codifferential matrix of $\delta$ 
using the definition of boundary and expressing the right-hand side $\delta\bar y$ in the basis of the $n$-simplices.

\heading{Persistence of the primary obstruction---the algorithm.} 
We remind that in the persistent setting the input contains
a~filtration of simplicial complexes $X\supseteq A_{r_0} \supseteq A_{r_1},\ldots,\supseteq A_{r_h}$ 
and a simplicial map $f': A_{r_0}\to \Sigma^{n-1}$.
We want to compute the smallest~value $j$ such that the restriction of $f'$ to ${A_{r_j}}$ can be extended to $A_{r_j}\cup X^{(n)}$.

Let $A:=A_{r_0}$ and $z,y$ be defined as above.
The cochain extension $\bar y$ of $(f')^\sharp(z)$ is also an extension of $(f'|_{A'})^\sharp(z)$ for each $A'\subseteq A$.
Thus we fix one $\bar y$ for all spaces $A_r$.
Then the only thing that is changing in solving $\delta c=\delta \bar y$, $c\in C^{n-1}(X,A_r;\Z)$ with increasing $r$,
is that $X\setminus A_r$ is becoming larger and we are allowed to include more columns into our matrix representing $\delta$.

Now we describe the algorithm on a lower level.
First we choose the cocycle $z\in\Sigma^{n-1}$ that generates the $(n-1)$ cohomology and its pullback
$y:=(f'|_A)^\sharp(z)\in Z^{n-1}(A;\Z)$ for $A=A_{r_0}$. 
%This is a cochain that assigns
%$1$ (or $-1$) to simplices $[v_1,\ldots,v_n]$ such that $(f'(v_1),\ldots,f'(v_n))$ is an even (or odd)
%permutation of $(e_1,\ldots,e_n)$, and $0$ if $\{f'(v_1),\ldots,f'(v_n)\}\neq\{e_1,\ldots,e_n\}$.
\item We fix an arbitrary extension $\bar y\in C^{n-1}(X;\Z)$ of $y$: the simplest option is to choose 
$\bar y(\Delta)=0$ for all $(n-1)$-simplices $\Delta\in X$ that are not in $A$.
\item We compute the filtration values of all $(n-1)$ simplices in $X$. 
%The filtration $r(\Delta)$ of $\Delta$ is the largest $r$
%so that $\Delta\in A_r$. If the filtration $A_{r_j}$ is induced by (\ref{e:filtration}), then equation 
by (\ref{e:filtration}). 
%directly defines the filtration $r(\Delta)$ of $\Delta$.
\item We order the $(n-1)$-simplices of $X$ by their filtration value, and choose an arbitrary enumeration of the $n$-simplices.
These choices will serve as bases of the $(n-1)$- and $n$-cochains (we identify simplices and their characteristic cochains).
\item We construct the matrix $M$ representing the codifferential with respect to the bases chosen above. The columns of $M$ are 
coboundaries of the $(n-1)$-simplices ordered by increasing filtration values.
%(Equivalently, we may compute $M$ as the transpose of the boundary matrix, which is computationally simpler.)
Further, we convert the right-hand side $\delta \bar y$ to an integer~vector $\aa$ using the chosen basis of $n$-simplices.
Recall that we want to solve $\delta c=\delta\bar y$ for $c$ that is a linear combination of $(n-1)$-simplices 
with filtration values smaller than $r_j$ for $j$ as small as possible. This directly translates to the following problem, which is the last
step of the persistence-of-primary-obstruction algorithm.

{\bf Problem {\sc Earliest Solution}}
Input: A matrix $M\in \Z^{p\times q}$ and a column vector $\aa\in \Z^p$.
Output: A column vector $\xx\in \Z^q$ such that $M\xx = \aa$.
Objective: Minimize the index of the \emph{lowest nonzero entry} of $\xx$, that is, $\ell\ge 0$ such that $x_\ell\neq 0$ and $x_{\ell+1}= x_{\ell+2} =\ldots = x_q = 0$.}

This problem could be solved by binary search on the value $\ell$ while solving an~ordinary 
linear system of equations in each iteration. 
Our implementation uses a simple matrix reduction approach (resembling algorithms for persistent homology) 
which avoids the binary search (see Appendix~\ref{s:pers} for details).
%\footnote{It is provably better in the case of
%equations over a~finite field.  In the case of integers, a coefficient blowup could make the running time exponential in the worst case. 
%However, we noticed no  significant growth of coefficients in our instances where the matrix $M$ has a rather special structure.} 
%The implementation details of \textsc{Earliest Solution} are in Appendix~\ref{s:pers}.

\heading{Secondary obstruction---extendability to $A\cup X^{(n+1)}$.}
Computing the secondary obstruction and its persistence contains similar ingredients but is more technical and we postpone a lower-level
description to Appendix~\ref{s:secondary-persistence}. Here we outline the main steps for the non-persistent version with a fixed $A$.
We assume that $f': A\to\Sigma^{n-1}$ is extendable to $A\cup X^{(n)}$ and that $\Omega(A)$ (described by (\ref{e:E})) is nonempty.

We need to implement the ``Steenrod square'' operation $\smile_{n-3}: C^{n-1}(X;\Z_2)\times C^{n-1}(X;\Z_2)\to C^{n+1}(X;\Z_2)$ 
on the level of cocycles. 
Peter FRANEK's avatar
Peter FRANEK committed
The algorithm for $\smile_{n-3}$ directly follows from formulas in~\cite[p. 292--293]{Steenrod}. 
Peter FRANEK's avatar
Peter FRANEK committed
674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708
For the following facts, we refer to \cite{Steenrod} and \cite{Steenrod:CohomologyOperationsObstructionsExtendingContinuousFunctions-1972}:
\begin{proposition}[Secondary obstruction]
Let $A\subseteq X$ be a pair of simplicial complexes, $f': A\to\Sigma^{n-1}$ be simplicial and 
assume that the ordering of vertices of $X$ and $\Sigma^{n-1}$ is chosen so that
$v\preceq w\,\,\Rightarrow\,\,f'(v)\preceq f'(w)$.
For each $x\in\Omega(A)$ let $(x\mod 2)$ be the image of $x$ under the natural homomorphism $C^{n-1}(X;\Z)\to C^{n-1}(X;\Z_2)$.  
v(x):=(x\mod 2)\smile_{n-3} (x\mod 2)
vanishes on $A$, that is, it is an~element of $Z^{n+1}(X,A;\Z_2)$.

Further, if $n>3$, then $f'$ can be extended to a map $X^{(n+1)}\to \Sigma^{n-1}$ iff $v(x)$ 
is a~relative coboundary for some $x\in\Omega(A)$.
Thus extendability to $X^{(n+1)}$ is equivalent to satisfiability of the equation $\delta c=v(x)$, $c\in C^n(X,A;\Z_2)$, 
for \emph{some} $x\in\Omega(A)$. To decide this, we parametrize $\Omega(A)$ by a fixed representant $x$ and generators $g_j$ of 
$Z^{n-1}(X,A;\Z)$: arbitrary element of $\Omega(A)$ is then $x-\sum_j u_j g_j$ for some $u_j\in\Z$.
To reduce the number of $j$'s, we only need to take generators of the cohomology group $H^{n-1}(X,A;\Z)$.
Exploiting the ``linearity'' of the operation $v$  on the level of cohomology (\cite[p. 504]{Steenrod}), we have that 
$v(x-\sum_j u_j g_j)$ is a coboundary iff $v(x)-\sum_j u_j v(g_j)$ is a coboundary. Thus our equation reduces to
$\delta c + \sum_j u_j v(g_j)=v(x)$.
This is a~system of equations with right-hand side $v(x)$ and unknowns $c$ and $u_j$, this time over the $\Z_2$-coefficients.

We also remark that the last proposition is valid also in the case $n=3$ once we replace $\Z_2$-coefficients by $\Z$-coefficients,
but deciding satisfiability of $\delta c=v(x)$ for some $x$ is hard (and undecidable for general spaces $X$).
%because in this dimension, the operation $\smile_{n-3}$ does not induce a~$\Z$-linear map on the level of cohomology,
%but rather a quadratic one. 
We show at the end of Appendix~\ref{s:secondary-persistence} 
that if $X$ is a triangulation of the topological cube $[0,1]^4$ and $n=3$, then
triviality of the secondary obstruction can easily be tested as well: this case is also included in our implementation.

Peter FRANEK's avatar
Peter FRANEK committed
709 710 711
To compute the smallest $r$ such that the map $f'|_{A_r}$ is extendable to $A_r\cup X^{(n+1)}$, we could use a binary search.
As in the case of the primary obstruction, it can be avoided and we can compute the persistence of the secondary obstruction
using a~single matrix reduction: this is explained in Appendix~\ref{s:secondary-persistence}.
Peter FRANEK's avatar
Peter FRANEK committed
712 713 714

\section{Application for robust optimization}
Peter FRANEK's avatar
Peter FRANEK committed
715 716 717 718
One application of zero verification is in optimization problems in which the solution of
$f(x) = 0$ is the approximation of the domain on which we want to maximize a~given objective function. 
Usually, the uncertainty has been modeled by unknown parameters $p$ in
the constraint $f(x, p)=0$, or $f(x,p)\leq 0$. These parameters are taken from a given set
Peter FRANEK's avatar
Peter FRANEK committed
719 720
(an ellipsoid or polygon) and a~common goal is to approximate the \emph{robust counterpart}, 
a~worst-case maximum of the objective function~\cite{Ben:2009}.
Peter FRANEK's avatar
Peter FRANEK committed
721 722 723

Our obstruction algorithm can be used if $f$ is uncertain and we don't have access to any finitely parametrized space of candidates.
Let $u: X\to\R$ be an objective function and assume that we want to approximate $\max\, u(x)$ subject to the constraint
Peter FRANEK's avatar
Peter FRANEK committed
724 725 726 727 728 729
$f(x)=0$.\footnote{In many optimization problems, the domain is defined by \emph{equations and inequalities} $f(x)=0\,\wedge \,g(x)\geq 0$. 
However, we can to some extent get rid of the inequalities $g_j(x)\geq 0$ by restricting the domain to the subcomplex $X'\subseteq X$ on which all components 
$g_j$'s are provably positive, resp. all $r$-perturbations of $g_j$ are provably positive.} 
As before, $f$ is given via function values on vertices of $X$ and a~simplexwise Lipschitz constant 

Peter FRANEK's avatar
Peter FRANEK committed
730 731
A number $\beta$ is a lower bound for the maximum of $u$ if and only if $f$ has a~zero on $u^{-1}[\beta,\infty)$. 
In general, $u^{-1}[\beta,\infty)$ is not a simplicial subcomplex
Peter FRANEK's avatar
Peter FRANEK committed
of $X$, so it is natural to discretize it.
Peter FRANEK's avatar
Peter FRANEK committed
733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750
For any $\beta\in \R$, we define $u_\beta$ to be the simplicial subcomplex 
of $X$ containing all simplices $\sigma$ such that $u(x)\geq \beta$ for all $x\in \sigma$. 

The filtration $(u_\beta)_{\beta\in\R}$ is step-like, $u_\beta\subseteq u^{-1}[\beta,\infty)$ and
we assume that $u_\beta$ is computable.\footnote{It is clearly computable in the important cases when $u$ is simplexwise linear.}
Then we can find a maximal $\beta$ for which we can guarantee a~zero of $f$ on
$u_\beta$: such $\beta$ is a lower bound on $\max_{x\in f^{-1}(0)} u(x)$.

On a~lower level, we can choose $r > \alpha n^{1/p}$ so that, by Lemma~\ref{t:approx}, $f'|_{A_{r}}$ is homotopic to $f|_{A_{r}}$.
Using the extendability oracle, we can decide, for each $\beta$, whether $f'|_{A_{r}\cap u_\beta}$ is extendable to all of $u_\beta$.
Via a~binary search\footnote{In fact, the binary search can be avoided within the primary/secundary obstruction computation.
For instance, the primary obstruction is measured by solvability of the equation $\delta c=\delta\bar{y}$ where 
$c\in Z^{n-1}(u_{\beta},A_r\cap u_{\beta})$.
If we include, in the coboundary matrix, only columns corresponding to $(n-1)$-simplices with filtration value $< r$ and
order the rows ($n$-simplices) by the $u$-filtration, and then
perform a column matrix reduction such that, after the reduction, the lowest nonzero element in each column is the last nonzero element in that row,
then the desired $\beta$ is the filtration value corresponding to the row of the lowest nonzero element on the right hand side after the reduction.},
we compute the largest $\beta$ such that $f'|_{A_{r}\cap u_\beta}$ is not extendable: by Section~\ref{sec:algorithm-oracle}
Peter FRANEK's avatar
Peter FRANEK committed
the robustness of zero of $f$ on $u_{\beta}$ is at least $r-\alpha>0$ and 
Peter FRANEK's avatar
Peter FRANEK committed
\beta \leq \inf_{\|g-f\|\leq r-\alpha} \,\,\max_{x\in g^{-1}(0)} u(x)\leq \max_{x\in f^{-1}(0)} \, u(x).
Peter FRANEK's avatar
Peter FRANEK committed
756 757 758 759 760 761

More generally, we could approximate the function
M(r):=\inf_{\|g-f\|\leq r} \,\max_{x\in g^{-1}(0)} u(x)
as follows. 
Peter FRANEK's avatar
Peter FRANEK committed
762 763
Define $U_\beta$ to be the simplicial subcomplex of $X$
containing all faces of all simplices $\sigma$ such that $u(x)\geq \beta$ for \emph{some} $x\in \sigma$.
Peter FRANEK's avatar
Peter FRANEK committed
We have the inclusions 
Peter FRANEK's avatar
Peter FRANEK committed
765 766 767 768 769
$$u_\beta\subseteq u^{-1}[\beta,\infty)\subseteq U_\beta.$$
For a given $\beta$, we may use our extendability oracle to compute a number $r(\beta)$ such that $f$ has an $r(\beta)$-robust zero on 
$u_\beta$ and a number $R(\beta)$ such that
$f$ has no $R(\beta)$-robust zero on $U_\beta$. 
The robustness of zero on $u^{-1}[\beta,\infty)$ is between $r(\beta)$ and $R(\beta)$ and
Peter FRANEK's avatar
Peter FRANEK committed
Peter FRANEK's avatar
Peter FRANEK committed
\min \{ \beta: \,\, r\in [r(\beta), R(\beta)] \}
Peter FRANEK's avatar
Peter FRANEK committed
772 773 774 775 776 777 778
Peter FRANEK's avatar
Peter FRANEK committed
\max \{\beta: \,\, r\in [r(\beta), R(\beta)] \}.
Peter FRANEK's avatar
Peter FRANEK committed
780 781 782 783

\section{Experimental results}\label{s:experiment}
Peter FRANEK's avatar
Peter FRANEK committed
Peter FRANEK's avatar
Peter FRANEK committed
785 786 787 788 789
One motivation for implementing the algorithm was to experimentally analyse the following question: 
\item How typical is a~situation in which the zero cannot be detected by primary obstruction and higher obstructions are needed?

Peter FRANEK's avatar
Peter FRANEK committed
790 791 792 793
To illustrate the flavour of this problem, consider a~function $f$ from an $(n+1)$-ball $B^{n+1}$ to $\R^n$ such that
the zero set is a~circle. If $r$ is small, then the $r$-neighborhood of the zero set is homeomorphic to a~solid torus $S^1\times B^n$.
An $n$-hyperplane intersecting the zero set transversally will typically 
intersect this torus in a $n$-disc $\{*\}\times B^n$ with a~zero of $f$ inside: this reflects the non-extendability to the $n$-skeleton.
Peter FRANEK's avatar
Peter FRANEK committed
However, with increasing $r$ (and hence increasing our freedom to perturbate the function), 
the primary obstruction will die once the $r$-neighborhood will touch the boundary or become a full $(n+1)$-ball:
Peter FRANEK's avatar
Peter FRANEK committed
796 797 798 799 800 801 802 803
in the latter case, a~nontrivial secondary obstruction is reflected by the homotopy class of the map from the boundary of this 
$(n+1)$-ball to $S^{n-1}$. This homotopy class is encoded in the gradient-induced framing of the original zero set of $f$: if the framing 
is trivial (framed null-cobordant), then higher obstruction don't occur. If the framing is ``twisted'', then they do.

Intuitively, we assumed that using Gaussian random fields, the gradient-induced framing of the zero set should be quite random and we would observe
twisted as well as untwisted cases. Experiments, however, do not support this so far, which we find surprising.

\heading{Description of the computation experiments.}
Peter FRANEK's avatar
Peter FRANEK committed
The lowest-dimensional case where nontrivial secondary obstruction can occur is $\dim X=4$ and $n=3$. 
Peter FRANEK's avatar
Peter FRANEK committed
805 806
Using an experimental approach, we generated random continuous functions from a~regular $4$-dimensional cubical grid into $\R^3$ 
taken from different probability distributions. 
Peter FRANEK's avatar
Peter FRANEK committed
807 808
The space $X$ was either a $4$-cube or a $4$-torus $(S^1)^4$ and
the underlying simplicial complex was the Freudenthal triangulation of the canonical cubical subdivision of $X$~\cite[p. 154]{allgower2003introduction}.
Peter FRANEK's avatar
Peter FRANEK committed
Instead of $A_r$ from Definition~\ref{d:A-r}, we used a courser filtration $A_r^\square$ based on the cubical structure,
Peter FRANEK's avatar
Peter FRANEK committed
see Appendix~\ref{s:implementation} for details.
Peter FRANEK's avatar
Peter FRANEK committed
We computed the vertex-approximation $f'$ from Definition~\ref{d:A-r} and the smallest $r_0>0$ such that
Peter FRANEK's avatar
Peter FRANEK committed
812 813
$f'$ is simplicial on $A_{r_0}^\square$. Then we found the persistence of the primary obstruction $r_1\geq r_0$ and the persistence
of the secondary obstruction $r_2\geq r_1$: the goal was to check whether instances with $r_2>r_1$ occur and how often. 
Peter FRANEK's avatar
Peter FRANEK committed

Peter FRANEK's avatar
Peter FRANEK committed
815 816 817 818 819
First we experimented with Gaussian random fields. 
Such functions are continuous and infinitely differentiable~\cite[Sec. 2.2]{Adler:1981}. 
For each component $f_i$ of $f$ and each vertex $x$, 
the random variable $f_i(x)$ was normalized to the standard normal distribution $N(0,1)$ 
and the covariance between $f_i(x)$ and $f_i(y)$ was taken to be $C(x,y)=\tilde{C}(|x-y|)$: we tried different functions $C$.
820 821
First we generated random functions such that the 
discrete Fourier transform of $C(0,x)$ was proportional to $((1+|p|^2)^{-l})_{p\in \{0,\ldots, g-1\}^4}$ for various constants $l$ 
Peter FRANEK's avatar
Peter FRANEK committed
(compare~\cite[p. 12]{Lang:2011}). The value $l=0$ corresponds to white noise and $l=\infty$ to constant functions.
Peter FRANEK's avatar
Peter FRANEK committed
823 824
While this procedure naturally creates functions on a torus, for experiments on a cube  we generated a~random function
on the discrete torus $\{1,2,\ldots, 2g\}^4$ and restricted it to the coordinates $\{1,\ldots, g\}^4$ to avoid periodicity.
Peter FRANEK's avatar
Peter FRANEK committed
825 826 827
The three components of $f$ were generated independently. To assure that the resulting function has zero at all, we analyzed the function
$f(x)-f(x_0)$ instead of $f(x)$, where $x_0$ was the midpoint of the cube, resp. a~fixed point in the torus.

Peter FRANEK's avatar
Peter FRANEK committed
In most cases, we detected a nontrivial primary obstruction, but not a single instance with secondary obstruction $r_2>r_1$.
Peter FRANEK's avatar
Peter FRANEK committed
829 830
To give an illustration, the following table shows some statistics of one of the experiments on a $4$-cube:
$l$ is the parameter of the distribution,
Peter FRANEK's avatar
Peter FRANEK committed
831 832 833 834
$g$ is the number of vertices in each dimension, 
$r_0$ the smallest value for which $f'|_{A_{r_0}^\square}$ is simplicial,
$r_1$ the average persistence of the primary obstruction in cases when $r_1>r_0$, 
and max. $r_1$ the largest persistence of primary obstruction. 
Peter FRANEK's avatar
Peter FRANEK committed
The averages are taken out of 1000 functions for $l\in \{3, 3.5,4, 4.5\}$ and out of 10\,\,000 for $l=5.0$.
Peter FRANEK's avatar
Peter FRANEK committed
Peter FRANEK's avatar
Peter FRANEK committed

Peter FRANEK's avatar
Peter FRANEK committed
838 839 840 841 842 843 844 845
$l$ & $g$ & $r_0$ & \% of $r_1>r_0$ & average $r_1$ if nontrivial & max. $r_1$\\
3.0  & 30 & $0.4$ & 78\% & $0.56$ & $0.95$ \\
3.5  & 30 & $0.21$  & 91\% & $0.41$ & $0.82$\\
4.0  & 25 & $0.15$ & 91\% & $0.32$ & 0.66 \\
4.5  & 25 & $0.1$ & 91\% & $0.25$ & $0.55$ \\
5.0  & 20 & $0.1$ & 87\% & $0.21$ & $0.63$ \\
Peter FRANEK's avatar
Peter FRANEK committed
\end{tabular} \\
Peter FRANEK's avatar
Peter FRANEK committed

Peter FRANEK's avatar
Peter FRANEK committed
When performing such experiments on the $4$-torus, it sometimes happened that the cup square of a computed extension $x\in\Omega(A_{r_1})$
Peter FRANEK's avatar
Peter FRANEK committed
was nontrivial in $H^4(X,A_{r_1}^\square)$, giving some ``hope'' of a nontrivial secondary obstruction: however, in all cases, this could be removed 
Peter FRANEK's avatar
Peter FRANEK committed
850 851 852 853 854
after replacing $x$ by another extension of the pullback $y=(f'|_A)^*(z)$ to the $2$-skeleton 
(see Section~\ref{s:obstructions}).\footnote{In fact, nontriviality
of the secondary obstruction on a $4$-torus can only be reduced to a~system of quadratic Diophantine equations. While we cannot 
algorithmically check satisfiability of quadratic equations, 
in all cases where we had to deal with this problem, these equations were almost trivial and solvable.}
Peter FRANEK's avatar
Peter FRANEK committed

Peter FRANEK's avatar
Peter FRANEK committed
In other rounds of experiments, we generated functions from a $5$-torus into $\R^4$ or
Peter FRANEK's avatar
Peter FRANEK committed
replaced the correlation function $C(x,y)$ by the Gaussian function
Peter FRANEK's avatar
Peter FRANEK committed
$\exp\left(-\frac{|x-y|^2}{2l^2}\right)$ for suitable $l>0$, but the results were were similar to that from the distribution above.
Peter FRANEK's avatar
Peter FRANEK committed
859 860 861 862 863

One possible explanation for the lack of secondary obstruction is that the cohomology in dimension $4$ has typically 
lower persistence than in dimension $3$ and most generators have already died when the primary obstruction (element of~$H^3$)
dies. A similar phenomenon occurs in persistent homology of excursion sets of random scalar fields, where the persistence
barcodes in dimension $0$ die before the barcodes in dimension 1, compare~\cite{Adler:2010}. 
Peter FRANEK's avatar
Peter FRANEK committed
In the vast majority of our experiments on the cube, the $4$-dimensional cohomology group $H^4(X,A_{r_1}^\square)$ 
Peter FRANEK's avatar
Peter FRANEK committed
865 866 867 868
was trivial for $r_1$ being the persistence of primary obstruction.
The lack of top dimensional cohomology in this case probably reflects the fact that most components of the neighborhood of the zero set intersect 
the boundary of the domain, although this argument does not apply for the torus. 

In another attempt to detect secondary obstruction in random fields we generated random homogenous quadratic polynomials on $[-1,1]^4$.
Peter FRANEK's avatar
Peter FRANEK committed
The coefficients $a_{ij}^k$ in $f_k(x)=\sum_{i,j} a_{i,j}^k x_i x_j$ were independent samples from a standard normal 
Peter FRANEK's avatar
Peter FRANEK committed
871 872 873 874 875 876 877
distribution.\footnote{This is motivated by the fact that the simplest examples of functions with nontrivial secondary obstruction are quadratic
and homogenous.}
The zero set of homogenous quadratic functions is either the origin alone or a cone intersecting the boundary $\partial [-1,1]^4$:
only the first case can yield a nontrivial $H^4(X,A_r^\square)$ and a nontrivial secondary obstruction. 
We generated around 70 thousand instances of random quadratic functions on a $10^4$ grid: around 2.2\% of them
had only the origin as the zero set, but there was no nontrivial secondary obstruction in a single instance. 

Peter FRANEK's avatar
Peter FRANEK committed
878 879 880 881

\heading{Experiments with formulas.}
Another motivation for implementing the algorithm was to test the running time and memory limitations in practice. 
Peter FRANEK's avatar
Peter FRANEK committed
882 883
Our testing benchmarks consisted of inputs in which the function values $f(v)$ were generated via formulas with known properties 
in a cubical grid. We ran many testing examples, some of them being shown in Appendix~\ref{s:formulas}.
Peter FRANEK's avatar
Peter FRANEK committed
884 885 886 887
To summarize the results, the performance is much better than the worst-case complexity bound derived in Section~\ref{s:complexity}
and is approximately linear in the number of simplices of the input. We were able to run benchmarks up to $\dim X=8$ for small grids,
such as $5^8$: the largest coboundary matrix for which we computed a nontrivial obstruction had 40 million columns.

Peter FRANEK's avatar
Peter FRANEK committed
888 889 890 891 892
%Another observation is the following. 
%If the function is generated by a ``nice'' (such as quadratic) formula and the grid is so coarse that 
%our theoretical analysis (based on Lemma~\ref{t:approx}) fails to certify a zero, the
%algorithm usually computes a nontrivial persistence of obstructions $r$ that is often a~surprisingly good 
%estimate of the robustness of zero of the underlying function.
Peter FRANEK's avatar
Peter FRANEK committed
893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910

The input size (and hence computational complexity) depends heavily on the encoding of the simplicial complex.
For example, we may 
specify the set of \emph{all} simplices, or the set of all top-dimensional simplices.\footnote{In other situations,
the input specifying the simplicial complex could be even smaller. 
One example is specifying the vertex set in $\R^m$ and assuming the Delaunay triangulation.}
Therefore we study \emph{parameterized complexity} as a function of the following parameters.
Let $m$ be the dimension of $X$ and $n\leq m$ the dimension of the target space $\R^n$.
We define $N$ to be the maximum of the number of $k$-simplices
for $k\in \{n-1,n,n+1\}$ and $V$ the number of vertices. 
In addition to specifying $X$, the input contains the function values in all vertices, that is, $n\times V$ numbers.

We present complexity bounds as a~function of $N,V,m, n$. 

\heading{Primary obstruction---complexity.}
Peter FRANEK's avatar
Peter FRANEK committed
We assume that the function values $f(v)$ at the vertices are all rational vectors and that we can compare their absolute values 
Peter FRANEK's avatar
Peter FRANEK committed
912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953
$|f(v)|$, $|f(w)|$ in unit time (these numbers may be roots of rational numbers for $\ell_p$ norms).
Then computing the vertex approximation $f'(v)$ for each vertex $v$ via~(\ref{d:vertex-app}) amounts to $O(nV)$ operations.
Computing the filtration of all $(n-1)$-simplices via formula (\ref{e:filtration}), as well as the pullback $y$ and its codifferential
 are by definition subroutines of complexity $O(nN)$; ordering the $(n-1)$-simplices by filtration is done in $O(N\log N)$.
The computation of the codifferential matrix is again of order $O(nN)$ if we store it in a~sparse format, 
because each row of the matrix corresponds to the boundary of an $n$-simplex and has only $n+1$ nonzero elements. 

The bottleneck of computing the primary obstruction is the EARLIEST SOLUTION algorithm described on page~\pageref{p:es}. 
An implementation based on a~binary or exponential search requires at most $\log N$ 
solutions of a linear system of Diophantine equations. Each of them is a system of at most $N$ rows and columns, all coefficients being $\pm 1$ or $0$.
By~\cite[Thm. 19]{storjohann1996fast+} we may solve such Diophantine system in $O(N^4\,\log^4 N)$ time, which yields 
$O(N^4\,\log^5 N+n(N+V))$ as a complexity bound for the primary obstruction. 
Assuming the lack of blowup of matrix coefficients during the matrix reduction, we can bound the number of arithmetic operations in EARLIEST SOLUTION  by $O(N^3)$. This is discussed in more detail in Appendix~\ref{s:pers}. 
%In practice, coefficients grow sublinear for examples of $N$ up to ten millions that we have tried.
In this scenario, sub-cubic bounds could be achieved using randomization~\cite[Thm. 39]{storjohann2005shifted}.
In practice, however, our implementation
of~EARLIEST SOLUTION exhibits subquadratic scaling, allowing us to experiment with instances
for $N \le 10^7$. This is not entirely surprising---large instances of simplicial boundary matrices
are commonly reduced in the field of computational topology. 

\heading{Secondary obstruction---complexity.}
The bottleneck in the secondary obstruction algorithm is the computation of all Steenrod squares of all generators of $H^{n-1}(X,A_r;\Z)$ 
for all filtration values $r$. In a naive implementation we may compute a set of generators of $Z^{n-1}(X;\Z)$ 
and their respective filtration values.
Generators of the kernel (over $\Z$) of a~matrix with at most $N$ rows and columns can be computed in 
$O(N^4)$~\cite[Theorem 1]{Buchmann:99}.
The number of such generators is bounded by $N$. In the Steenrod square computation,
we need to compute, in the worst case, the values on all $(n+1)$-simplices; in each evaluation, the formula for $\smile_{n-3}$
described in~\cite{Steenrod} contains an iteration over all elements of $m\choose 4$ (Steenrod pairs).
Thus, computing the Steenrod squares of the generators of $Z^{n-1}(X;\Z)$ is $O(N^2\,m^4)$.
The final matrix computation corresponding to equation (\ref{e:second-pers3}) is done over the field $\Z_2$ 
which only requires $O(N^\omega)$ operations for a constant $\omega<3$~\cite[Proposition 6]{jeannerod2013rank}.
This yields a~complexity bound of $O(N^4+N^2 m^4)$ for the persistence of secondary obstruction. For all practical purposes,
it is safe to assume that the values of $m$ can be neglected.

%\bibliographystyle{spbasic}      % basic style, author-year citations
\bibliographystyle{spmpsci}      % mathematics and physical sciences
%\bibliographystyle{spphys}       % APS-like style for physics
%\bibliography{}   % name your BibTeX data base

Peter FRANEK's avatar
Peter FRANEK committed
Peter FRANEK's avatar
Peter FRANEK committed
955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135 1136 1137 1138 1139 1140 1141 1142 1143 1144 1145 1146 1147 1148 1149 1150 1151 1152 1153 1154 1155 1156 1157 1158 1159 1160 1161 1162 1163 1164 1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216 1217 1218 1219 1220 1221 1222 1223 1224 1225 1226 1227 1228 1229 1230 1231 1232 1233 1234 1235 1236 1237 1238 1239 1240 1241 1242 1243 1244 1245 1246 1247 1248 1249 1250 1251 1252 1253 1254 1255 1256 1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269 1270 1271 1272 1273 1274 1275 1276 1277 1278 1279 1280 1281 1282 1283 1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297 1298 1299 1300 1301 1302 1303 1304 1305 1306 1307 1308 1309 1310 1311 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324 1325 1326 1327 1328 1329 1330 1331 1332 1333 1334 1335 1336 1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399 1400 1401 1402 1403 1404 1405 1406 1407 1408 1409 1410 1411 1412 1413 1414 1415 1416 1417 1418 1419 1420 1421 1422 1423 1424 1425 1426 1427 1428 1429 1430 1431 1432 1433 1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444 1445 1446 1447 1448 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461 1462 1463 1464 1465 1466 1467 1468 1469 1470 1471 1472 1473 1474 1475 1476 1477 1478 1479 1480 1481 1482 1483 1484


\section{Secondary obstruction.}
\heading{Persistence of the secondary obstruction---the algorithm for $n>3$.} 
Assume that a filtration $X\supseteq A_{r_0}\supseteq A_{r_1}\supseteq \ldots$ and a simplicial map $f': A_{r_0}\to\Sigma^{n-1}$
are given, $n>3$, and vertices on $X$ and $\Sigma^{n-1}$ are ordered so that $f'$ is order-preserving: this order is used
in the implementation of the $\smile_{n-3}$ operation on the level of cochains.
Further, we assume that the persistence of primary obstruction $r_j$ has already been computed by the algorithm
described on page~\pageref{page:primary-persistence}.
%That is, the restriction of $f$ to $A_{r_j}$ is extendable to some continuous map $A_{r_j}\cup X^{(n)}\to\Sigma^{n-1}$. 
We continue to use the notation of Section~\ref{s:obstructions}: in particular, $z$ is the characteristic cocycle of a~fixed
$(n-1)$-simplex in $\Sigma^{n-1}$, $\bar y\in C^{n-1}(X,A_{r_j})$ is a cochain extending the pullback $y=(f')^\sharp(z)\in Z^{n-1}(A_{r_j})$ 
of $z$ and 
$\emptyset\neq \Omega(A_{r_j})$ is the set of all $(n-1)$-cocycles on $X$ that extend $y$ on $A_{r_j}$. 

By Proposition~\ref{p:secondary}, the persistence of secondary obstruction is the smallest number $r_k$ such that
\delta c=v(x) \quad \text{for some } c\in C^n(X,A_{r_k};\Z_2), x\in \Omega(A_{r_k})
\end{equation} where $v$ is defined by (\ref{e:steenrod}).

Let $x\in\Omega(A_{r_j})$ be a~fixed extension of $\bar y$, computed in the algorithm for primary persistence. 
Then also $x\in\Omega(A_{r_k})$ for each $k\geq j$.
For any $k$, $\Omega(A_{r_k})$ is a coset in $Z^{n-1}(X)$
and hence equation (\ref{e:second-pers1}) reduces to
\delta c=v(x-w), \quad \text{for some } w\in Z^{n-1}(X,A_{r_k};\Z), \,\, c\in C^n(X,A_{r_k};\Z_2).
The crucial property that we will use is that $v(x-w)$ is a relative coboundary iff $v(x)-v(w)$ is a relative coboundary:
this follows directly from the linearity of the Steenrod square operation $H^{n-1}(X,A;\Z_2)\to H^{n+1}(X,A;\Z_2)$ for $n>3$.
Thus we can reformulate (\ref{e:second-pers2}) to the problem of finding the minimal $r_k$ such that
\delta c+v(w) = v(x) \quad \text{for some } w\in Z^{n-1}(X,A_{r_k};\Z), \,\, c\in C^n(X,A_{r_k};\Z_2).
To simplify the computations, we don't need to consider all cocycles $w\in Z^{n-1}(X,A_r)$ 
but only generators of the cohomology $H^{n-1}(X,A_r;\Z)$: the Steenrod square of any relative coboundary is again a relative 
coboundary $\delta c'$ for some $c'\in C^{n}(X,A_r;\Z_2)$, so adding it has no impact on the solvability of (\ref{e:second-pers3}).

The right-hand side of (\ref{e:second-pers3}) $v(x)$ is a cocycle that does not 
depend on $k$. The left-hand side is a~combination of coboundaries of characteristic cocycles of $n$-simplices $\Delta^n$
and cochains of the form $v(w)$ for $(n-1)$-cocycles $w$. To each $\Delta$ and $w$ is assigned a filtration value and 
we want to minimize the maximal filtration $r_k$ of all the $\Delta^n$'s and $w$'s.

Summarizing the above steps, we obtain the following algorithm:
\item Order the vertices  of $\Sigma^{n-1}$ and the vertices of $X$ so that $f$ is order-preserving
\item For a precomputed $x\in\Omega(A_{r_j})$, compute the relative cocycle $v(x)\in Z^{n+1}(X,A_{r_j};\Z_2)$ by the 
definition of the Steenrod operation $\smile_{n-3}$
\item Compute a subset $W\subseteq Z^{n-1}(X)$ that contains all cohomology generators $w$ of all $H^{n-1}(X, A_r;\Z)$ for all $r\geq r_j$.
(It may be a~set of generators of $Z^{n-1}(X)$.) To each $w\in W$, assign its filtration value $r(w)$ 
to be the minimal $r$ such that $w$ is zero on $A_r$.
\item Compute the filtration value of all $n$-simplices, using (\ref{e:filtration}). Order 
\item Order all $n$-simplices and elements of $W$ by their filtration.
\item Choose a basis of $(n+1)$-simplices and express the right-hand side $v(x)$ as a~vector $\aa\in (\Z_2)^q$. 
\item Create a matrix $M$ whose column set consists of
\item all coboundaries of characteristic cochains of $n$-simplices
expressed in the basis of $(n+1)$-simplices over $\Z_2$
\item all elements $v(w)$ expressed in the basis of $(n+1)$-simplices
\item Order the columns of $M$ by filtration.
Computing the persistence of the secondary obstruction then reduces to solving the EARLIEST SOLUTION problem for $M\xx = \aa$, 
this time over the $\Z_2$-field.

The hardest part is to compute the cohomology generators: this algorithm is summarized as follows.

{\bf Problem {\sc Persistent Generators:}}
Input: A filtration $X\supseteq A_0\supseteq\ldots\supseteq A_h=\emptyset$.
Output: A sequence $g_1,\ldots,g_\nu\in Z^{n-1}(X;\Z)$ and a sequence of integers $\mu(0)\le\ldots\le\mu(h)$ such that $[g_1],\ldots, [g_{\mu(i)}]$ generate $H^{n-1}(X,A_i;\Z)$ for each $i$.
%Objective: Minimize the index of the \emph{lowest nonzero entry}
%of $\xx$, that is, $\ell\ge 0$ such that $x_\ell\neq 0$ and $x_{\ell+1}=
%x_{\ell+2} =\ldots = x_q = 0$.
We give imeplementation details of this part on a lower-level in Section~\ref{s:persistent-gen}.
Our algorithm will give the output with the number of generators $\nu$ minimal in a certain sense. Notably, when working over $\Q$ or $\Z_p$ instead of integers, the output would correspond to a persistence bar code with all the death information erased but including representative (co)cycles for each bar.
\heading{The special case $f: [0,1]^4\to\R^3$.}
This section justifies the claim that for $X=[0,1]^4$ and $n=3$,
the general algorithm for computing persistence of secondary obstruction works, if we completely ignore the persistence generators $w$
and perform all computations over $\Z$-coefficients.
The following text is not completely self-contained, as the proof of Lemma~\ref{l:trivial-cup} requires some definitions and techniques from
the field of algebraic topology. 

Assume first that $n=3$ and $X$ is arbitrary. The two main differences compared to the algorithm above are:
\item The Steenrod operation and final coboundary matrix have to be computed over $\Z$, not $\Z_2$. This goes back to the fact that
the homotopy group $\pi_3(S^2)\simeq \Z$, unlike $\pi_n(S^{n-1})\simeq \Z_2$ for $n>3$. The homotopy group serves as cohomology 
coefficients in the theory of obstructions.
\item The operation $\smile_{n-3}$ reduces to the cup product $\smile$ and the operation $w\mapsto w\smile w$ is not 
linear on the level of cohomology but quadratic.
This means that while for any particular extension $x\in \Omega(A)$ of $\bar y$, we may test satisfiability of $\delta c=v(x)$,
we cannot test the \emph{existence} of such $x$ using a linear system of equations.
However, if $X=[0,1]^4$, we claim that for any two extensions $x,y\in\Omega(A)$ of $\bar y$, $v(x)-v(y)$ is a relative coboundary
and thus we need to check the equation $\delta c=v(x)$ only for one $x$.
Let $X$ be a triangulation of $[0,1]^4$, $A\subseteq X$, $x\in Z^{2}(X)$ and $w\in Z^2(X,A)$. Then 
(x-w)\smile (x-w) \,-\, (x\smile x) \in B^4(X,A).
Using the parametrization $\Omega(A)=\{x-w:\,\,w\in Z^{n-1}(X,A;\Z)\}$ we immediately obtain that our general algorithm
for the persistence of secondary obstruction works once we replace the $\Z_2$-coefficients by $\Z$-coefficients in its final step. 
Moreover, we may completely ignore the persistent generators $w$ and don't need to compute them at all.\footnote{Note that
$x\smile x$ may represent a nontrivial element of $H^4(X,A)$, as $x$ is not an element of $Z^2(X,A)$ in general. The fact that
$x\smile x$ is zero on $A$ follows from the fact that $f$ is order-preserving.}
By bilinearity of the cup product, $(x-w)\smile (x-w) \, - \, x\smile x\,= \, - \, x\smile w \, -\, w\smile x\, +\, w\smile w$.
The mixed-term $x\smile w$ is a relative coboundary, because $\smile$ induces a bilinear product on the level of cohomology
$H^2(X)\times H^2(X,A)\to H^4(X,A)$ and $H^2(X)=0$, as $X$ is contractible. It remains to show that $w\smile w$ is a coboundary.
Let $CA$ be the cone over $A$ and $A\hookrightarrow CA$ be an inclusion.
The inclusion $A\hookrightarrow X$ can be extended to a map $CA\to X$, because $X$ is contractible. The map of pairs
$(CA,A)\to (X,A)$ induces the commutative diagram in which the rows are the long exact sequences of cohomology groups.
H^{*-1}(X) & \to & H^{*-1}(A) & \to & H^*(X,A) &\to & H^*(X) &\to & H^*(A) \\
\downarrow & &\downarrow & &\downarrow& &\downarrow& & \downarrow\\
H^{*-1}(CA) & \to & H^{*-1}(A) & \to & H^*(CA,A) &\to & H^*(CA) &\to & H^*(A) 
The vertical arrows $H^{*(-1)}(X)\to H^{*(-1)}(CA)$ are trivial as both spaces are contractible and 
$H^{*(-1)}(A)\to H^{*(-1)}(A)$ are identities.
By the five-lemma~\cite[p. 129]{Hatcher}, the middle homomorphism $H^*(X,A)\to H^*(CA,A)$ is an isomorphism. 
Further, $H^j(CA,A)\simeq H^j(CA/A)$ for $j>0$. The space $CA/A=:\Sigma A$ is the \emph{suspension} of $A$.
The cup product $H^2(\Sigma A)\times H^2(\Sigma A)\to H^{4}(\Sigma A)$ 
of a~suspension is trivial~\cite[Corollary 4.11]{Bredon-at} and the naturality of cup product implies that the cup product 
$H^2(X,A)\times H^2(X,A)\to H^{4}(X,A)$ is trivial for $j>0$ as well. This shows that $w\smile w\in B^4(X,A)$. 

\section{Persistent integral homology computations}\label{s:pers}
\subsection{Algorithm for the \textsc{Earliest Solution} problem}

The {\tt earliest\_solution} algorithm is used to find the persistence of a (co)cycle. 
This is closely related to computing persistent homology, which is a well-studied problem, at least for coefficient in a finite field. 
We adapt the boundary matrix reduction algorithm by 
Edelsbrunner, Letscher and Zomorodian~\cite{edelsbrunner2002}, originally developed for persistent homology. Note that this algorithm, unlike classical 
Gaussian elimination or Smith normal form algorithms, is incremental, which is required by our application. 
Moreover, it only uses column operations, making efficient implementation for sparse matrices relatively easy.

\heading{Efficiency over finite fields.}
Recent work on computing persistent homology over finite fields resulted in significant performance improvements, see~\cite{phat,gudhi}. Despite the cubic worst-case bound, linear scaling is achieved on practical datasets, involving sparse matrices of size $10^9 \times 10^9$ and more. This encouraged us to adapt the modern version of the 
classical persistence algorithm to solve our problem over the integers, rather than adapting 
classical Smith or Hermite normal form algorithms to the persistent setting.

\heading{Reduced form and reduction.}
We adapt the notation common in computational topology literature: the \emph{lowest nonzero} of a nonzero column is defined as the lowest position (largest index) with nonzero coefficient. A sub-matrix is called \emph{reduced} if all the lowest nonzeros are unique. In other words, given a lowest nonzero, there may be other nonzero entries in the same row, but they must not be \emph{lowest} nonzeros. When this invariant is not satisfied, we say there is a collision. By \emph{lowest value}, we refer to the value of the lowest nonzero entry.

The algorithm starts from an empty matrix and adds one column at a time, maintaining the reduced prefix of the matrix, $R$. The rightmost column of each prefix is called the current column. 

Procedure {\tt reduce\_column} reduces the column {\tt curr} with respect to the reduced prefix {\tt R}.

def reduce_column(current, R, force_divisibility=False):        
        while curr collides with some column in R:
                coll = column in R colliding with curr
                P = lowest_value(coll)
                Q = lowest_value(curr)             
                if force_divisibility is True and P does not divide Q:
                        return curr
                use Euclid to find nonzero a,b,c,d s.t. 
                    gcd(a,b) == 1 and gcd(c,d) == 1 and
                    a*P + b*Q == gcd(P,Q) and
                    c*P + d*Q == 0                          
                coll, curr = a * coll + b * curr, c * coll + d * curr                
        return curr

Procedure {\tt earliest\_solution} solves the stated problem.
def earliest_solution(M, a):
        augment M with identity matrix (above)
        augment a with zeros (above)
        R = []
        for col in M:
                reduced = reduce_column(col, R, force_divisibility=False)
                append reduced to R            
                reduce_column(a, R, force_divisibility=True)            
                if a is zero:
                        return change of basis of a 
        return None

The basic operation is the addition of two different columns, without affecting the column span of the relevant matrix prefix. Therefore the solution is unaffected.

To retrieve the solution, at each step we attempt to reduce the input column vector {\tt a} with respect to the currently reduced prefix {\tt R}. The solution exists iff {\tt a} becomes zero, and is encoded by (the negation of) the change of basis column of {\tt a}. 
Since we solve the equation $Ax = a$ (and not $Ax = ka$) we perform additional divisibility check:
see {\tt force\_divisibility} in the {\tt reduce\_column} procedure.

\heading{From finite fields to integers.}
In the integral case, we may modify both the current column and the colliding column. This is in contrast with the finite field case, in which 
only the current column is modified. To determine the required linear combination of columns, we use the extended Euclid algorithm. 
%This is clearly more efficient than adding the entire columns multiple times to solve a single collision. 
As a result, the lowest \emph{value} of a certain column might change (decrease) many times during the reduction of other columns, but the \emph{position} is fixed once the column is reduced. Because of this, while reducing column {\tt a}, we need to take into account 
previously reduced columns, and not only the current one. Moreover, after a colliding column
is affected, it may not be in the column span of the prefix of the original matrix ending at this column. At this stage, however, this column may only affect columns that succeed the currently reduced column. Therefore, correctnes of the algorithm is unaffected.

%In the algorithm (for finite fields) the current column is reduced by adding multiples of columns from the prefix. More precisely, whenever the current column causes a collision, a multiple of the colliding column is added, so that the lowest nonzero of the column is higher (has smaller index), or the column is empty (all coefficients are zero). The key property is that the lowest nonzeros are unique (while the reduced matrix depends on reduction scheme). 

\heading{Efficiency.} For efficiency reasons we use one technique suggested in~\cite{phat}: The current column is stored in a data-structure handling fast column additions and maximum element queries. One natural choice is a balanced binary search tree; 
more efficient alternatives are available. This way we avoid the following common bad case: Let $n$ be the total number of columns in $M$. When an current column becomes dense (the number nonzero entries is $\Theta(n)$), adding $\Theta(n)$ sparse columns takes time quadratic in $n$. Avoiding this situation does \emph{not} imply that we can efficiently handle matrices that become dense due to fill-in. However in practice, often a small number of columns display this behavior.

We also perform the computations in an on-line fashion -- we read columns one by one and stop once a solution is found. This gives significant practical improvements, because often the necessary matrix prefix is very small compared to the matrix of the entire complex.

Overall the algorithm performs well, exhibiting roughly linear scaling in the number of nonzero entries of the original matrix. In particular, we didn't observe coefficient blowup. 
Note that the lowest nonzero position in the current column decreases after resolving each collision, so the number of collision resolutions is quadratic in $n$. Each such resolution requires combining two columns, potentially of size $O(n)$, due to fill-in. This leads to cubic worst case running time bound, assuming that the magnitude of coefficients can be bounded by a constant. Of course, there exist cases when the blowup does occur, and the above analysis does not apply. More advanced algorithms can be used to alleviate the effect of the blowup, possibly at the cost of simplicity and efficiency in the situations that we encountered thus far. 

\subsection{Algorithm for the \textsc{Persistent Generators}
After possibly refining the sequence $(A_i)_i$ we may assume that for
each filtration value $i$ there is exactly one $(n-1)$-simplex
$\Delta_i^{n-1}$ of that filtration value, that is $\Delta_i^{n-1}
\in A_{i-1}\setminus A_{i}$ (and possibly several $(n-2)$-simplices
in $A_{i-1}\setminus A_{i}$). This makes the algorithm and its analysis simpler.

The algorithm for \textsc{Persistent Generators} problem follows.
By the statement ``reduce the column'' we refer to the procedure
\texttt{reduce\_column} from Section~\ref{s:earliest}
\item Let $M$ be the coboundary matrix for $\delta^{n-1}\:C^{n-1}(X;\Z)\to
C^n(X;\Z)$, that is, it consists of columns $\mm$ for integral
combination for $\delta \Delta^{n-1}$ for each $\Delta^{n-1}\in
X$ of growing filtration values. In exactly the same way we create
the coboundary matrix $N$ for $\delta^{n-2}$. During the reduction
of the matrix $M$ we keep track of the change of basis. (On low
level, we augment the matrix $M$ and then the change of basis
vector for a given column is encoded in the augmented part of
that column.)

\item We initialize $G:=()$ to be an empty sequence of $(n-1)$-cocycles
on $X$ and $\mu(-1):=0$.
\item For each filtration value $i=0,\ldots,h$ do the following:
\item Reduce all columns of $N$ of the filtration value $i$.
\item Reduce the column $\mm$ of $M$ of filtration value $i$
(i.e., corresponding to the cochain $\delta \Delta^{n-1}_i$).
Let $\gg$ be the change of basis vector after the reduction.
\item If the reduced column $\mm$ equals to $0$ and there is
no reduced column $\nn$ in $N$ such that $\lowest(\nn)=i$ and
$\nn_{\lowest(\nn)}|\gg_{\lowest(\gg)}$ then add $\gg$ to $G$
and set $\mu(i):=\mu(i-1)+1$.
\item Otherwise set $\mu(i):=\mu(i-1)$
\item Output $G$ and $\mu(1),\ldots,\mu(h)$.

The above algorithm solves the \textsc{Persistent Cycles} problem. Namely, after its $i$th iteration, the cohomology group $H^{n-1}(X,A_i)$ is generated by $[g_1],\ldots,[g_{\mu(i)}]$ where cycles $g_1,\ldots,g_{\mu(i)}\in Z^{n-1}(X,A_i;\Z)$ correspond to the columns
$\gg$ added to $G$ during the iterations $0,\ldots,i$.
We proceed by induction on the value $i$. For $i=-1$ the claim trivially holds.

Let us assume that $i\ge 0$ and that the claim holds for $i-1$. First observe that once the new column $\mm$ is not reduced to a zero column, then $Z^{n-1}(X,A_i;\Z)=Z^{n-1}(X,A_{i-1};\Z)$. Otherwise $Z^{n-1}(X,A_i;\Z)=Z^{n-1}(X,A_{i-1};\Z)\oplus \langle g\rangle$ where the cocycle $g\in Z^{n-1}(X,A_i;\Z)$ corresponds to the change of basis vector $\gg$ in the $i$th iteration. We only have to check whether $[g]$ is a linear combination of  $[g_1], \ldots, [g_{\mu(i-1)}]$. By the induction hypothesis, it happens if and only if $g$ is cohomologous to a cocycle in $Z^{n-1}(X,A_{i-1};\Z)$. And this in turn is true if and only if we can reduce to zero the lowest nonzero component (that is, the $i$th component) of the vector $\gg$ by adding a combination of columns from $N$ of filtration value at most $i$ (these columns generate the group of coboundaries $B^{n-1}(X,A_{i};\Z)$). Since this is exactly the reduced part of the matrix $N$, it is enough to find if there is a column with the lowest nonzero index equal to $i$ and check the divisibility condition.  
\marek{Alternatively, we could describe the algorithm by a pseudocode:}

def persistent_cycles(M,N):
    augment M
    RM = {}, RN = {}
    cycles = []
    for each filtration value i:        
        for each column col of  N with filtration value i: 
            reduced = reduce_column(col,RN)
            append reduced to RN
        col = the column column in M  with filtration value i
        reduced = reduce_column(col,RM)
        append reduced to RM
        if reduced is zero:
            g=augmented part of reduced
            if there is no column K in RN  s.t.
                        lowest(K)==lowest(g) and
                        lowest_value(K) divides lowest_value(g):
                append g to cycles
                m[i] = m[i-1] + 1
                m[i] = m[i-1]                
     return cycles

\section{Some details of our implementation: exploiting the cubical structure.}
The domain $X$ we consider in the implementation is a cube $[0,1]^m$ triangulated as follows. 
We define $I_1,\ldots, I_m$ to be unit intervals subdivided into $n_i$ equidistant intervals of length $1/n_i$, $i=1,\ldots, m$.
This yields a cubical set structure $\prod_j I_j$ on the unit cube. Further, we subdivide each $m$-cube into $m!$ simplices
via the Freudenthal triangulation~\cite[p. 154]{allgower2003introduction}.  The resulting triangulation
naturally corresponds to the product $\prod I_j$, understood as a product $\prod_j I_j$ of \emph{simplicial sets}~\cite{Friedm08}.
We call the set of vertices a \emph{grid}. 
A function $f:X\to\R^n$ is then given by a set of $\R^n$-vectors in each vertex (a~multidimensional rasterized image), 
together with a~simplexwise Lipschitz constant $\alpha$.

Many operations in the algorithm outlined in the paper body---such as the computation of the pullback $y=f^\sharp(z)$, 
computing codifferentials of cochains and the $\smile_{n-3}$ operation---can be done \emph{locally}, 
without having to work with the full lists of simplices of a given dimension. 
The only step where global structure is needed is the {\textsc{Earliest Solution}} algorithm, 
where the matrix $M$ has rows indexed by the set of all $k$-simplices and columns indexed by $(k-1)$-simplices.
Indeed, the construction of this matrix and computation with it are the most memory- and time-consuming operations,
as the number of $k$-simplices increases rapidly with dimension.

Therefore, we implemented a modification of the algorithm described above based on the more feasible \emph{cubical} structure.
Instead of working with the simplicial filtration $A_r$, we use the filtration $A_r^\square\subseteq A_s^\square$ where 
$A_r^\square$ is defined to be the triangulation of the set of all cubes $c$ such that $|f(v)|\geq v$ for all vertices $v$ of $c$. 
These are still simplicial complexes and $f'$, $y$, the extension $\tilde{y}$ of $y$ and $\delta\tilde{y}$ are defined with no changes. For the computation
of persistence, however, we switch to the cubical setting via the Eilenberg-Zilber reduction~\cite{Eilenberg:1953}. 
We denote by $C^*(\Pi_j I_j)$ the simplicial cochains and by $\otimes_j C^*(I_j)$ the cubical cochains: there exist
chain homomorphisms of degree $0$
&AW: \otimes_j C^*(I_j)\to C^*(\Pi_j I_j)\\
&EML: C^*(\Pi_j I_j)\to \otimes_j C^*(I_j)
such that $EML\circ AW$ is the identity and $AW\circ EML$ is chain homotopic to the identity. Both maps induce cohomology 
isomorphisms and can be relatively easily implemented using common formulas~\cite{Pedro:2005}. This allows us to switch between the simplicial
and cubical cochains anytime we need. Within the computation of persistence of the primary obstruction,
we compute the smallest $j$ so that there exists a cubical cochain $c_\square\in Z^{n-1}(X,A_{r_j}^\square)$ such that
$\delta c_\square=(\delta\tilde{y})_\square:=EML (\delta \tilde{y})$.
The matrix computation part \textsc{Earliest Solution} deals with the cubical coboundary
matrix, which is significantly smaller than the simplicial one.
For an illustration, the number of $k$ simplices in the triangulation of one $m$-cube is 50 already for $(m,k)=(4,2)$ and
more than 4 millions for $(m,k)=(8,6)$. 

For the secondary obstruction, we need to convert $c_\square$ back into a simplicial cochain and construct a simplicial $(n-1)$-cochain 
$x\in \Omega(A_{r_j})$ that extends $y$ on $A_{r_j}$ and is a global cocycle. 
We denote by SHI the cochain homotopy map $SHI: C^*(\Pi_j I_j)\to C^{*-1} (\Pi_j I_j)$, 
satisfying $AW\circ EML - \mathrm{id}=\delta\circ SHI + SHI\circ\delta$.
Then, using $\delta^2=0$, we obtain
\delta AW c_\square 
&= AW \delta c_\square = AW (\delta\tilde{y})_\square = AW\, EML\,\delta\tilde{y}= \\
& = \delta AW \,\,EML\, \tilde{y}=\delta \tilde{y} + \delta \,SHI \,\delta \tilde{y}.
We compute $\tilde{c}:=AW (c_\square) - SHI(\delta \tilde{y})$ which is a \emph{simplicial} cochain zero on $A_{r_j}^\square$
and the last equation asserts that $\delta \tilde{c}=\delta\tilde{y}$ which allows us to compute the simplicial cocycle
x:=\tilde{y} - \tilde{c} \in \Omega(r_j)
having avoided to work with large lists of all simplices. 
The computation of $v(x):=(x\mod 2)\smile_{n-3} (x\mod 2)$ is done on the simplicial level and the property 
$x\in Z^{n+1}(X,A_{r_j}^\square ;\Z_2)$ depends on the fact that the chosen ordering of vertices of $X$ and $\Sigma^{n-1}$
is compatible with $f'$. Then we again apply the $EML$ operator to $v(x)$ and convert it to a cubical cochain $v(x)_\square$.
Persistent generators $w\in Z^{n-1}(X,A_{r}^\square;\Z)$ and their Steenrod images can be computed on the cubical level
and the cubical persistence of the secondary obstruction is done via a~cubical coboundary matrix, this time in dimension $n,(n+1)$.

There is a price to pay for the more convenient cubical filtration: it is courser than the simplicial one and
we can no more use the estimate on the robustness of zero set 
derived in Section~\ref{sec:algorithm-oracle}. On one hand, we have the relation $A_r^\square\subseteq A_r$.
This implies that whenever $f'|_{A_r}$ is homotopic to $f/|f|$, then so is $f'|_{A_r^\square}$, and
non-extendability of $f'|_{A_r^\square}$ to higher skeleta of $X$ implies the non-extendability of $f'|_{A_r}$.
Thus whenever the algorithm certifies non-extendability of $f|_{A_r^\square}$ and $r>\alpha n^{1/p}$, then 
$r-\alpha$ is a lower bound on the robustness of zero.
If we can additionally prove that $A_s\subseteq A_r^\square$, then the extendability of $f'|_{A_r^\square}$
implies the extendability of $f'|_{A_s}$.  The simplexwise $\alpha$-Lipschitz condition implies that any two vertices
$u,v$ of a cube $c$ satisfy $|f(u)-f(v)|\leq 2\alpha$, as $u,v$ can be connected by  path of at most two simplicial edges.
This implies $A_{r+2\alpha}\subseteq A_r^\square$ and extendability of $f'|_{A_r^{\square}}$ to all of $X$ implies the upper bound 
$r+3\alpha$ on the robustness.

Based on the experience with various functions in low dimensions, the practical performance has not one significant bottleneck.
The most resources-consuming steps in the primary obstruction computation include the computation of the cubical filtration and
the EARLIER SOLUTION subroutine, especially if the matrix has millions of columns. 

\section{Experiments with formulas.}
Our prototype is implementation in Python using numpy. While the efficiency is severely limited by this choice, 
we were able to run examples for $\dim X\leq 8$.
%Even in these cases, the estimates derived from Lemma~\ref{t:approx} are often so limiting that the smallest $r$
%satisfying them is already larger than the robustness. 

Although the strength of our methods is primarily in cases where we have uncertainty on $f$, 
the following examples provide some observations about the performance in cases where the
function values at the vertices are computed using exact formulas. 

In what follows, we assume a subdivision of the interval $[-1,1]$ into a set $I_g$ of $g$ equidistant points
and consider a grid $I_g^m\subseteq [-1,1]^m$ of size $g^m$. We chose the max-norm $\ell_\infty$ on $\R^n$ which yields
the smallest value of $\alpha n^{1/p}$ and allows smaller ``initial $r_0$''. This is desired, as for small grids it often happens
that $\alpha$ is large and $\alpha n^{1/p}$ may be larger then the robustness of zero in which case we fail to detect anything.

\heading{Testing primary obstruction on a quadratic function.}
First we used the function $f: [-1,1]^n\to\R^n$ given by 
& f_1(x)=x_1^2-x_2^2-\ldots -x_n^2 \\
& f_2(x)=2x_1 x_2\\
& \ldots \\
& f_n(x)=2x_1 x_n.
This function has a single zero in the origin and has been used as a benchmark example in~\cite{nondec}. 
If $n$ is even, then the zero is robust and its robustness equals 
$\min_{x\in\partial [-1,1]^m} |f(x)|$: the zero has index two in this case. If $n$ is odd, then the zero has index zero 
and can be removed by arbitrary small perturbations. 
Using the max-norm, we use calculus to derive the estimate 
$|f(x)-f(y)|_{\infty}\leq 2n\,|x-y|_\infty$. Whenever $x,y$ are contained in one simplex of the triangulation, then the max-norm satisfies
$|x-y|_\infty\leq 2/(g-1)$ and we can use the piecewise Lipschitz constant $\alpha:=4n/(g-1)$. 

For our experiments, we now assume that $f$ is only given via a~list of function values in the grid $g^n$ and the 
simplexwise Lipschitz constant $\alpha=4n/(g-1)$. In this case, we have $\dim X=n$ so only primary obstruction can be nontrivial.
Using Lemma~\ref{t:approx}
for $p=\infty$ yields that whenever $r>\alpha$, then $f'$ is simplicial on $A_r$ and homotopic to $f/|f|$.
The relation $A_r^\square\subseteq A_r$ then implies that $f'$ is simplicial and homotopic to $f/|f|$ on $A_r^\square$ as well.
Thus we can run our algorithm with an initial $r_0=\alpha=4n/(g-1)$, compute the persistence of primary obstruction $r_1\geq r_0$
and---using the estimates at the end of Section~\ref{s:implementation}---conclude that whenever $r_1>r_0$ 
then the robustness of zero is between $r_1-\alpha$ and $r_1+3\alpha$.

The following table illustrates some properties of the computations: $n$ refers to the dimension, $g$ is the 
number of points subdividing the intervals in each dimension, the initial $r_0=\alpha$ is as above, 
$r_1$ is the persistence of the primary obstruction, 
$\delta^{n-1}$ the number of columns of the cubical coboundary matrix $C^{n-1}_\square\to C^n_\square$, 
the ``computed robustness'' displays lower and upper bounds $r_1-\alpha$, $r_1+3\alpha$ on the robustness of zero 
and the ``true robustness'' column the approximation of the real robustness of zero of the function $f$ 
defined exactly by the above formula,\footnote{This is $0$ in odd dimensions and $\frac{2}{\sqrt{n}+1}$ in even dimension.}
all wrt. the $|\cdot |_\infty$ norm.
$n$ & $g$ & $r_0:=\alpha$ & $r_1$ & \# columns of $\delta^{n-1}$ & computed robustness & true robustness\\
$2$ & $20$ & $0.421$ & $0.8643$ & $760$ & $[0.443,2.127]$ & $0.828$\\
$2$ & $100$ & $0.08$ & $0.8285$ & $19 \,\,800$ & $[0.748,1.07]$ & $0.828$\\
$2$ & $500$ & $0.016$ & $0.83$ & $499\,\,000$ & $[0.814,0.878]$ & $0.828$\\
$3$ & $20$ & $0.632$ & $=r_0$ & $21\,\,660$ & $\leq 2.379$ & $0$\\
$3$ & $50$ & $0.245$ & $=r_0$ & $360\,\,150$ & $\leq 0.971$ & $0$\\
$3$ & $100$ & $0.121$ & $=r_0$ & $2\,\,940\,\,300$ & $\leq 0.483$ & $0$\\
$4$ & $20$ & $0.842$ & $=r_0$ & $548\,\,720$ & $\leq 3.37$ & $0.667$\\
$4$ & $30$ & $0.552$ & $0.711$ & $2\,\, 926\,\, 680$ & $[0.159, 2.367]$ & $0.667$\\
$4$ & $40$ & $0.41$  & $0.667$ & $9\,\,491\,\,040$  &  $[0.257, 1.897]$ & $0.667$
The total running times of these 9 computations is displayed in Figure~\ref{fig:times}.
\caption{Total running time of the computations, performed on Intel Xeon CPU X5680, 3.3 GHz,
as a function of the number of $(n-1)$-simplices (which equals the number of columns in the coboundary matrix).
%\vskip -2cm
The data suggest that the at least in these simple cases, the running time is approximately linear in the number of $(n-1)$-simplices.
The $n=3$ case is below the other two, because the primary obstruction is trivial there and the 
EARLIEST SOLUTION subroutine terminates almost immediately, using only very few columns in the matrix reduction.

In higher dimensions, the condition $r_0\geq \alpha:=4n(g-1)$ is too strict and we can only hardly start with $r_0$ that is smaller than the robustness.
However, we could take $r_0$ to be the minimal value for which $f'|_{A_{r_0}^\square}$ is simplicial:
in all cases that we have tried, we verified that $f$ was nowhere zero on $A_{r_0}^\square$ with $f'$ homotopic to $f/|f|$ 
on $A_{r_0}^\square$. 

Surprisingly, starting with minimal $r_0$ for which $f'|_{A_{r_0}^\square}$ is simplicial and computing 
the minimal $r_1$ for which $f'$ is extendable, yields in many cases a~much better estimate of the robustness of zero set than
the estimates based on Lemma~\ref{t:approx}: this can already be seen in the above table where the $r_1$ column is a good approximation
of the true robustness. 
We give the table with smaller grids that continues to higher dimensions.
$n$ & $g$ & min simplicial $r_0$ & $r_1$ & true robustness\\
$2$ & $10$ & $0.025$ & $0.889$ & $0.828$\\
$3$ & $10$ & $0.025$ & $=r_0$ & $0$\\
$4$ & $10$ & $0.025$ & $0.667$ & $0.667$\\
$5$ & $10$ & $0.074$ & $=r_0$ & $0$ \\
$6$ & $10$ & $0.074$ & $0.667$ & $0.58$ \\
$7$ & $6$ & $0.241$ & $=r_0$ & $0$ \\
$8$ & $5$ & $0.251$ & $1.0$ & $0.522$

It is an interesting question to find natural conditions on functions, other than the relatively weak simplexwise Lipschitz property,
that justify the usage of smaller grids and guarantee that the robustness of zero set is close to the computed minimal $r$ 
for which $f'|_{A_{r}^\square}$ becomes extendable. 

\heading{A function with nontrivial secondary obstruction.}
A second function we experimented with is $h: [-1,1]^{n+1}\to\R^n$ given by
& h_1(x)= 2x_0 x_2 + 2x_1 x_3\\
& h_2(x)= 2x_1 x_2 - 2x_0 x_3\\
& h_3(x)= x_0^2+x_1^2-x_2^2-x_3^2\\
& h_4(x)=x_4 \\
& \hskip 1cm \ldots \\
& h_n(x)=x_n.
The restriction of $h$ to $\partial [-1,1]^{n+1}$ is the generator of the nontrivial homotopy group $\pi_n(S^{n-1})$:
in case $n=3$ it is the Hopf map and for $n>3$ its iterated suspension.
The robustness of zero equals $\min_{x\in \partial [-1,1]^{n+1}} |h(x)|$ and\footnote{The robustness 
is $1$ in the $\ell_2$ norm and $\sqrt{3}-1$ in the max-norm.} 
it is the simplest example where a nontrivial secondary obstruction occurs. Common tests for zero verification such as
the degree test would fail here.

Again, we work with the max-norm $\ell_\infty$ and derive, via elementary calculus, the estimate $2(n+1)$ on the global Lipschitz constant.
This yields $|h(x)-h(y)|\leq \frac{4(n+1)}{g-1}=:\alpha$.
Assume that only $\alpha$ and a list of function values in a grid $g^{n+1}$ is given. If $g$ is large enough so that $r_0:=\alpha$
is smaller than $\min_{\partial} |h(x)|_{\infty}=\sqrt{3}-1$, then the algorithm computes the minimal $r_1>r_0$ for which
$h|_{A_r^\square}$ is extendable to a nowhere zero function $[-1,1]^{n+1}\to\R^{n}\setminus\{0\}$. For this to succeed, we need to take
$g$ at least 22 for $n=3$.

Again, the program gives surprisingly good results for much smaller grids, although Lemma~\ref{t:approx} gives no guarantees. 
In all cases that we tried, we verified that whenever $r_0$ is large enough so that $h'|_{A_{r}^\square}$ is simplicial, 
then it is homotopic to $h/|h|$ and the algorithm computes that the secondary obstruction dies close to the real robustness of zero.
$n$ & $g$ & min simplicial $r_0$ & $r_1$ & true robustness\\
$3$ & $10$ & $0.1$ & $0.79$ & $0.732$\\
$4$ & $10$ & $0.111$ & $0.79$ & $0.732$\\
$5$ & $8$ & $0.163$ & $0.816$ & $0.732$\\

% end of file template.tex