The complex Helmholtz equation $(\Delta + k^2)u=f$ (where $k\in{\mathbb R},u(\cdot),f(\cdot)\in{\mathbb C}$) is a mainstay of computational wave simulation. Despite its apparent simplicity, efficient numerical methods are challenging to design and, in some applications, regarded as an open problem. Two sources of difficulty are the large number of degrees of freedom and the indefiniteness of the matrices arising after discretisation. Seeking to meet them within the novel framework of probabilistic domain decomposition, we set out to rewrite the Helmholtz equation into a form amenable to the Feynman-Kac formula for elliptic boundary value problems. We consider two typical scenarios, the scattering of a plane wave and the propagation inside a cavity, and recast them as a sequence of Poisson equations. By means of stochastic arguments, we find a sufficient and simulatable condition for the convergence of the iterations. Upon discretisation a necessary condition for convergence can be derived by adding up the iterates using the harmonic series for the matrix inverse -- we illustrate the procedure in the case of finite differences. From a practical point of view, our results are ultimately of limited scope. Nonetheless, this unexpected -- even paradoxical -- new direction of attack on the Helmholtz equation proposed by this work offers a fresh perspective on this classical and difficult problem. Our results show that there indeed exists a predictable range $k<k_{max}$ in which this new ansatz works with $k_{max}$ being far below the challenging situation.
Binary codes of length $n$ may be viewed as subsets of vertices of the Boolean hypercube $\{0,1\}^n$. The ability of a linear error-correcting code to recover erasures is connected to influences of particular monotone Boolean functions. These functions provide insight into the role that particular coordinates play in a code's erasure repair capability. In this paper, we consider directly the influences of coordinates of a code. We describe a family of codes, called codes with minimum disjoint support, for which all influences may be determined. As a consequence, we find influences of repetition codes and certain distinct weight codes. Computing influences is typically circumvented by appealing to the transitivity of the automorphism group of the code. Some of the codes considered here fail to meet the transitivity conditions requires for these standard approaches, yet we can compute them directly.
Recently, Letzter proved that any graph of order $n$ contains a collection $\mathcal{P}$ of $O(n\log^\star n)$ paths with the following property: for all distinct edges $e$ and $f$ there exists a path in $\mathcal{P}$ which contains $e$ but not $f$. We improve this upper bound to $19 n$, thus answering a question of G.O.H. Katona and confirming a conjecture independently posed by Balogh, Csaba, Martin, and Pluh\'ar and by Falgas-Ravry, Kittipassorn, Kor\'andi, Letzter, and Narayanan. Our proof is elementary and self-contained.
We propose MNPCA, a novel non-linear generalization of (2D)$^2${PCA}, a classical linear method for the simultaneous dimension reduction of both rows and columns of a set of matrix-valued data. MNPCA is based on optimizing over separate non-linear mappings on the left and right singular spaces of the observations, essentially amounting to the decoupling of the two sides of the matrices. We develop a comprehensive theoretical framework for MNPCA by viewing it as an eigenproblem in reproducing kernel Hilbert spaces. We study the resulting estimators on both population and sample levels, deriving their convergence rates and formulating a coordinate representation to allow the method to be used in practice. Simulations and a real data example demonstrate MNPCA's good performance over its competitors.
Consider the problem of binary hypothesis testing. Given $Z$ coming from either $\mathbb P^{\otimes m}$ or $\mathbb Q^{\otimes m}$, to decide between the two with small probability of error it is sufficient and in most cases necessary to have $m \asymp 1/\epsilon^2$, where $\epsilon$ measures the separation between $\mathbb P$ and $\mathbb Q$ in total variation ($\mathsf{TV}$). Achieving this, however, requires complete knowledge of the distributions and can be done, for example, using the Neyman-Pearson test. In this paper we consider a variation of the problem, which we call likelihood-free (or simulation-based) hypothesis testing, where access to $\mathbb P$ and $\mathbb Q$ is given through $n$ iid observations from each. In the case when $\mathbb P,\mathbb Q$ are assumed to belong to a non-parametric family $\mathcal P$, we demonstrate the existence of a fundamental trade-off between $n$ and $m$ given by $nm \asymp n^2_\mathsf{GoF}(\epsilon,\cal P)$, where $n_\mathsf{GoF}$ is the minimax sample complexity of testing between the hypotheses $H_0: \mathbb P= \mathbb Q$ vs $H_1: \mathsf{TV}(\mathbb P,\mathbb Q) \ge \epsilon$. We show this for three families of distributions: $\beta$-smooth densities supported on $[0,1]^d$, the Gaussian sequence model over a Sobolev ellipsoid, and the collection of distributions on alphabet $[k]=\{1,2,\dots,k\}$ with pmfs bounded by $c/k$ for fixed $c$. For the larger family of all distributions on $[k]$ we obtain a more complicated trade-off that exhibits a phase-transition. The test that we propose, based on the $L^2$-distance statistic of Ingster, simultaneously achieves all points on the trade-off curve for the regular classes. This demonstrates the possibility of testing without fully estimating the distributions, provided $m\gg1/\epsilon^2$.
We consider the problem $(\rm P)$ of exactly fitting an ellipsoid (centered at $0$) to $n$ standard Gaussian random vectors in $\mathbb{R}^d$, as $n, d \to \infty$ with $n / d^2 \to \alpha > 0$. This problem is conjectured to undergo a sharp transition: with high probability, $(\rm P)$ has a solution if $\alpha < 1/4$, while $(\rm P)$ has no solutions if $\alpha > 1/4$. So far, only a trivial bound $\alpha > 1/2$ is known to imply the absence of solutions, while the sharpest results on the positive side assume $\alpha \leq \eta$ (for $\eta > 0$ a small constant) to prove that $(\rm P)$ is solvable. In this work we study universality between this problem and a so-called "Gaussian equivalent", for which the same transition can be rigorously analyzed. Our main results are twofold. On the positive side, we prove that if $\alpha < 1/4$, there exist an ellipsoid fitting all the points up to a small error, and that the lengths of its principal axes are bounded above and below. On the other hand, for $\alpha > 1/4$, we show that achieving small fitting error is not possible if the length of the ellipsoid's shortest axis does not approach $0$ as $d \to \infty$ (and in particular there does not exist any ellipsoid fit whose shortest axis length is bounded away from $0$ as $d \to \infty$). To the best of our knowledge, our work is the first rigorous result characterizing the expected phase transition in ellipsoid fitting at $\alpha = 1/4$. In a companion non-rigorous work, the first author and D. Kunisky give a general analysis of ellipsoid fitting using the replica method of statistical physics, which inspired the present work.
We consider the problem of approximating the solution to $A(\mu) x(\mu) = b$ for many different values of the parameter $\mu$. Here we assume $A(\mu)$ is large, sparse, and nonsingular with a nonlinear dependence on $\mu$. Our method is based on a companion linearization derived from an accurate Chebyshev interpolation of $A(\mu)$ on the interval $[-a,a]$, $a \in \mathbb{R}$. The solution to the linearization is approximated in a preconditioned BiCG setting for shifted systems, where the Krylov basis matrix is formed once. This process leads to a short-term recurrence method, where one execution of the algorithm produces the approximation to $x(\mu)$ for many different values of the parameter $\mu \in [-a,a]$ simultaneously. In particular, this work proposes one algorithm which applies a shift-and-invert preconditioner exactly as well as an algorithm which applies the preconditioner inexactly. The competitiveness of the algorithms are illustrated with large-scale problems arising from a finite element discretization of a Helmholtz equation with parameterized material coefficient. The software used in the simulations is publicly available online, and thus all our experiments are reproducible.
A widely used formulation for null hypotheses in the analysis of multivariate $d$-dimensional data is $\mathcal{H}_0: \boldsymbol{H} \boldsymbol{\theta} =\boldsymbol{y}$ with $\boldsymbol{H}$ $\in\mathbb{R}^{m\times d}$, $\boldsymbol{\theta}$ $\in \mathbb{R}^d$ and $\boldsymbol{y}\in\mathbb{R}^m$, where $m\leq d$. Here the unknown parameter vector $\boldsymbol{\theta}$ can, for example, be the expectation vector $\boldsymbol{\mu}$, a vector $\boldsymbol{\beta} $ containing regression coefficients or a quantile vector $\boldsymbol{q}$. Also, the vector of nonparametric relative effects $\boldsymbol{p}$ or an upper triangular vectorized covariance matrix $\textbf{v}$ are useful choices. However, even without multiplying the hypothesis with a scalar $\gamma\neq 0$, there is a multitude of possibilities to formulate the same null hypothesis with different hypothesis matrices $\boldsymbol{H}$ and corresponding vectors $\boldsymbol{y}$. Although it is a well-known fact that in case of $\boldsymbol{y}=\boldsymbol{0}$ there exists a unique projection matrix $\boldsymbol{P}$ with $\boldsymbol{H}\boldsymbol{\theta}=\boldsymbol{0}\Leftrightarrow \boldsymbol{P}\boldsymbol{\theta}=\boldsymbol{0}$, for $\boldsymbol{y}\neq \boldsymbol{0}$ such a projection matrix does not necessarily exist. Moreover, since such hypotheses are often investigated using a quadratic form as the test statistic, the corresponding projection matrices often contain zero rows; so, they are not even effective from a computational aspect. In this manuscript, we show that for the Wald-type-statistic (WTS), which is one of the most frequently used quadratic forms, the choice of the concrete hypothesis matrix does not affect the test decision. Moreover, some simulations are conducted to investigate the possible influence of the hypothesis matrix on the computation time.
Let $(X_t)_{t \ge 0}$ be the solution of the stochastic differential equation $$dX_t = b(X_t) dt+A dZ_t, \quad X_{0}=x,$$ where $b: \mathbb{R}^d \rightarrow \mathbb R^d$ is a Lipschitz function, $A \in \mathbb R^{d \times d}$ is a positive definite matrix, $(Z_t)_{t\geq 0}$ is a $d$-dimensional rotationally invariant $\alpha$-stable L\'evy process with $\alpha \in (1,2)$ and $x\in\mathbb{R}^{d}$. We use two Euler-Maruyama schemes with decreasing step sizes $\Gamma = (\gamma_n)_{n\in \mathbb{N}}$ to approximate the invariant measure of $(X_t)_{t \ge 0}$: one with i.i.d. $\alpha$-stable distributed random variables as its innovations and the other with i.i.d. Pareto distributed random variables as its innovations. We study the convergence rate of these two approximation schemes in the Wasserstein-1 distance. For the first scheme, when the function $b$ is Lipschitz and satisfies a certain dissipation condition, we show that the convergence rate is $\gamma^{1/\alpha}_n$. Under an additional assumption on the second order directional derivatives of $b$, this convergence rate can be improved to $\gamma^{1+\frac 1 {\alpha}-\frac{1}{\kappa}}_n$ for any $\kappa \in [1,\alpha)$. For the second scheme, when the function $b$ is twice continuously differentiable, we obtain a convergence rate of $\gamma^{\frac{2-\alpha}{\alpha}}_n$. We show that the rate $\gamma^{\frac{2-\alpha}{\alpha}}_n$ is optimal for the one dimensional stable Ornstein-Uhlenbeck process. Our theorems indicate that the recent remarkable result about the unadjusted Langevin algorithm with additive innovations can be extended to the SDEs driven by an $\alpha$-stable L\'evy process and the corresponding convergence rate has a similar behaviour. Compared with the previous result, we have relaxed the second order differentiability condition to the Lipschitz condition for the first scheme.
We expound on some known lower bounds of the quadratic Wasserstein distance between random vectors in $\mathbb{R}^n$ with an emphasis on affine transformations that have been used in manifold learning of data in Wasserstein space. In particular, we give concrete lower bounds for rotated copies of random vectors in $\mathbb{R}^2$ with uncorrelated components by computing the Bures metric between the covariance matrices. We also derive upper bounds for compositions of affine maps which yield a fruitful variety of diffeomorphisms applied to an initial data measure. We apply these bounds to various distributions including those lying on a 1-dimensional manifold in $\mathbb{R}^2$ and illustrate the quality of the bounds. Finally, we give a framework for mimicking handwritten digit or alphabet datasets that can be applied in a manifold learning framework.
We consider Gibbs distributions, which are families of probability distributions over a discrete space $\Omega$ with probability mass function of the form $\mu^\Omega_\beta(\omega) \propto e^{\beta H(\omega)}$ for $\beta$ in an interval $[\beta_{\min}, \beta_{\max}]$ and $H( \omega ) \in \{0 \} \cup [1, n]$. The partition function is the normalization factor $Z(\beta)=\sum_{\omega \in\Omega}e^{\beta H(\omega)}$. Two important parameters of these distributions are the log partition ratio $q = \log \tfrac{Z(\beta_{\max})}{Z(\beta_{\min})}$ and the counts $c_x = |H^{-1}(x)|$. These are correlated with system parameters in a number of physical applications and sampling algorithms. Our first main result is to estimate the counts $c_x$ using roughly $\tilde O( \frac{q}{\varepsilon^2})$ samples for general Gibbs distributions and $\tilde O( \frac{n^2}{\varepsilon^2} )$ samples for integer-valued distributions (ignoring some second-order terms and parameters), and we show this is optimal up to logarithmic factors. We illustrate with improved algorithms for counting connected subgraphs, independent sets, and perfect matchings. As a key subroutine, we also develop algorithms to compute the partition function $Z$ using $\tilde O(\frac{q}{\varepsilon^2})$ samples for general Gibbs distributions and using $\tilde O(\frac{n^2}{\varepsilon^2})$ samples for integer-valued distributions.