A widely used approach to compute the action $f(A)v$ of a matrix function $f(A)$ on a vector $v$ is to use a rational approximation $r$ for $f$ and compute $r(A)v$ instead. If $r$ is not computed adaptively as in rational Krylov methods, this is usually done using the partial fraction expansion of $r$ and solving linear systems with matrices $A- \tau I$ for the various poles $\tau$ of $r$. Here we investigate an alternative approach for the case that a continued fraction representation for the rational function is known rather than a partial fraction expansion. This is typically the case, for example, for Pad\'e approximations. From the continued fraction, we first construct a matrix pencil from which we then obtain what we call the CF-matrix (continued fraction matrix), a block tridiagonal matrix whose blocks consist of polynomials of $A$ with degree bounded by 1 for many continued fractions. We show that one can evaluate $r(A)v$ by solving a single linear system with the CF-matrix and present a number of first theoretical results as a basis for an analysis of future, specific solution methods for the large linear system. While the CF-matrix approach is of principal interest on its own as a new way to compute $f(A)v$, it can in particular be beneficial when a partial fraction expansion is not known beforehand and computing its parameters is ill-conditioned. We report some numerical experiments which show that with standard preconditioners we can achieve fast convergence in the iterative solution of the large linear system.
In the Hospitals/Residents problem, every hospital has an upper quota that limits the number of residents assigned to it. While, in some applications, each hospital also has a lower quota for the number of residents it receives. In this setting, a stable matching may not exist. Envy-freeness is introduced as a relaxation of stability that allows blocking pairs involving a resident and an empty position of a hospital. While, envy-free matching might not exist either when lower quotas are introduced. We consider the problem of finding a feasible matching that satisfies lower quotas and upper quotas and minimizes envy in terms of envy-pairs and envy-residents in the Hospitals/Resident problem with Lower Quota. We show that the problem is NP-hard with both envy measurement. We also give a simple exponential-time algorithm for the Minimum-Envy-Pair HRLQ problem.
We revisit the notion of root polynomials, thoroughly studied in [F. Dopico and V. Noferini, Root polynomials and their role in the theory of matrix polynomials, Linear Algebra Appl. 584:37--78, 2020] for general polynomial matrices, and show how they can efficiently be computed in the case of matrix pencils. The staircase algorithm implicitly computes so-called zero directions, as defined in [P. Van Dooren, Computation of zero directions of transfer functions, Proceedings IEEE 32nd CDC, 3132--3137, 1993]. However, zero directions generally do not provide the correct information on partial multiplicities and minimal indices. These indices are instead provided by two special cases of zero directions, namely, root polynomials and vectors of a minimal basis of the pencil. We show how to extract, starting from the block triangular pencil that the staircase algorithm computes, both a minimal basis and a maximal set of root polynomials in an efficient manner. Moreover, we argue that the accuracy of the computation of the root polynomials can be improved by making use of iterative refinement.
Let a polytope $\mathcal{P}$ be defined by one of the following ways: (i) $\mathcal{P} = \{x \in \mathbb{R}^n \colon A x \leq b\}$, where $A \in \mathbb{Z}^{(n+m) \times n}$, $b \in \mathbb{Z}^{(n+m)}$, and $rank(A) = n$, (ii) $\mathcal{P} = \{x \in \mathbb{R}_+^n \colon A x = b\}$, where $A \in \mathbb{Z}^{m \times n}$, $b \in \mathbb{Z}^{m}$, and $rank(A) = m$, and let all the rank minors of $A$ be bounded by $\Delta$ in the absolute values. We show that $|\mathcal{P} \cap \mathbb{Z}^n|$ can be computed with an algorithm, having the arithmetic complexity bound $$ O\bigl( \nu(d,m,\Delta) \cdot d^3 \cdot \Delta^4 \cdot \log(\Delta) \bigr), $$ where $d = \dim(\mathcal{P})$ and $\nu(d,m,\Delta)$ is the maximal possible number of vertices in a $d$-dimensional polytope $P$, defined by one of the systems above. Using the obtained result, we have the following arithmetical complexity bounds to compute $|P \cap \mathbb{Z}^n|$: 1) The bound $O(\frac{d}{m}+1)^m \cdot d^3 \cdot \Delta^4 \cdot \log(\Delta)$ that is polynomial on $d$ and $\Delta$, for any fixed $m$; 2) The bound $O\bigl(\frac{m}{d}+1\bigr)^{\frac{d}{2}} \cdot d^4 \cdot \Delta^4 \cdot \log(\Delta)$ that is polynomial on $m$ and $\Delta$, for any fixed $d$; 3) The bound $O(d)^{4 + \frac{d}{2}} \cdot \Delta^{4+d} \cdot \log(\Delta)$ that is polynomial on $\Delta$, for any fixed $d$. Given bounds can be used to obtain faster algorithms for the ILP feasibility problem, and for the problem to count integer points in a simplex or in an unbounded Subset-Sum polytope. Unbounded and parametric versions of the above problem are also considered.
We consider \emph{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 \emph{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 and perfect matchings in a graph. We develop a key subroutine to estimate the partition function $Z$. Specifically, it generates a data structure to estimate $Z(\beta)$ for \emph{all} values $\beta$, without further samples. Constructing the data structure requires $O(\frac{q \log n}{\varepsilon^2})$ samples for general Gibbs distributions and $O(\frac{n^2 \log n}{\varepsilon^2} + n \log q)$ samples for integer-valued distributions. This improves over a prior algorithm of Huber (2015) which computes a single point estimate $Z(\beta_\max)$ using $O( q \log n( \log q + \log \log n + \varepsilon^{-2}))$ samples. We show matching lower bounds, demonstrating that this complexity is optimal as a function of $n$ and $q$ up to logarithmic terms.
A general quantum circuit can be simulated classically in exponential time. If it has a planar layout, then a tensor-network contraction algorithm due to Markov and Shi has a runtime exponential in the square root of its size, or more generally exponential in the treewidth of the underlying graph. Separately, Gottesman and Knill showed that if all gates are restricted to be Clifford, then there is a polynomial time simulation. We combine these two ideas and show that treewidth and planarity can be exploited to improve Clifford circuit simulation. Our main result is a classical algorithm with runtime scaling asymptotically as $n^{\omega/2}<n^{1.19}$ which samples from the output distribution obtained by measuring all $n$ qubits of a planar graph state in given Pauli bases. Here $\omega$ is the matrix multiplication exponent. We also provide a classical algorithm with the same asymptotic runtime which samples from the output distribution of any constant-depth Clifford circuit in a planar geometry. Our work improves known classical algorithms with cubic runtime. A key ingredient is a mapping which, given a tree decomposition of some graph $G$, produces a Clifford circuit with a structure that mirrors the tree decomposition and which emulates measurement of the corresponding graph state. We provide a classical simulation of this circuit with the runtime stated above for planar graphs and otherwise $nt^{\omega-1}$ where $t$ is the width of the tree decomposition. Our algorithm incorporates two subroutines which may be of independent interest. The first is a matrix-multiplication-time version of the Gottesman-Knill simulation of multi-qubit measurement on stabilizer states. The second is a new classical algorithm for solving symmetric linear systems over $\mathbb{F}_2$ in a planar geometry, extending previous works which only applied to non-singular linear systems in the analogous setting.
In this paper we study the maximum degree of interaction which may emerge in distributed systems. It is assumed that a distributed system is represented by a graph of nodes interacting over edges. Each node has some amount of data. The intensity of interaction over an edge is proportional to the product of the amounts of data in each node at either end of the edge. The maximum sum of interactions over the edges is searched for. This model can be extended to other interacting entities. For bipartite graphs and odd-length cycles we prove that the greatest degree of interaction emerge when the whole data is concentrated in an arbitrary pair of neighbors. Equal partitioning of the load is shown to be optimum for complete graphs. Finally, we show that in general graphs for maximum interaction the data should be distributed equally between the nodes of the largest clique in the graph. We also present in this context a result of Motzkin and Straus from 1965 for the maximal interaction objective.
Given $n$ subspaces of a finite-dimensional vector space over a fixed finite field $\mathbb F$, we wish to find a "branch-decomposition" of these subspaces of width at most $k$ that is a subcubic tree $T$ with $n$ leaves mapped bijectively to the subspaces such that for every edge $e$ of $T$, the sum of subspaces associated to the leaves in one component of $T-e$ and the sum of subspaces associated to the leaves in the other component have the intersection of dimension at most $k$. This problem includes the problems of computing branch-width of $\mathbb F$-represented matroids, rank-width of graphs, branch-width of hypergraphs, and carving-width of graphs. We present a fixed-parameter algorithm to construct such a branch-decomposition of width at most $k$, if it exists, for input subspaces of a finite-dimensional vector space over $\mathbb F$. Our algorithm is analogous to the algorithm of Bodlaender and Kloks (1996) on tree-width of graphs. To extend their framework to branch-decompositions of vector spaces, we developed highly generic tools for branch-decompositions on vector spaces. The only known previous fixed-parameter algorithm for branch-width of $\mathbb F$-represented matroids was due to Hlin\v{e}n\'y and Oum (2008) that runs in time $O(n^3)$ where $n$ is the number of elements of the input $\mathbb F$-represented matroid. But their method is highly indirect. Their algorithm uses the nontrivial fact by Geelen et al. (2003) that the number of forbidden minors is finite and uses the algorithm of Hlin\v{e}n\'y (2006) on checking monadic second-order formulas on $\mathbb F$-represented matroids of small branch-width. Our result does not depend on such a fact and is completely self-contained, and yet matches their asymptotic running time for each fixed $k$.
In this work, we consider the distributed optimization of non-smooth convex functions using a network of computing units. We investigate this problem under two regularity assumptions: (1) the Lipschitz continuity of the global objective function, and (2) the Lipschitz continuity of local individual functions. Under the local regularity assumption, we provide the first optimal first-order decentralized algorithm called multi-step primal-dual (MSPD) and its corresponding optimal convergence rate. A notable aspect of this result is that, for non-smooth functions, while the dominant term of the error is in $O(1/\sqrt{t})$, the structure of the communication network only impacts a second-order term in $O(1/t)$, where $t$ is time. In other words, the error due to limits in communication resources decreases at a fast rate even in the case of non-strongly-convex objective functions. Under the global regularity assumption, we provide a simple yet efficient algorithm called distributed randomized smoothing (DRS) based on a local smoothing of the objective function, and show that DRS is within a $d^{1/4}$ multiplicative factor of the optimal convergence rate, where $d$ is the underlying dimension.
This paper describes a suite of algorithms for constructing low-rank approximations of an input matrix from a random linear image of the matrix, called a sketch. These methods can preserve structural properties of the input matrix, such as positive-semidefiniteness, and they can produce approximations with a user-specified rank. The algorithms are simple, accurate, numerically stable, and provably correct. Moreover, each method is accompanied by an informative error bound that allows users to select parameters a priori to achieve a given approximation quality. These claims are supported by numerical experiments with real and synthetic data.
In this paper, we study the optimal convergence rate for distributed convex optimization problems in networks. We model the communication restrictions imposed by the network as a set of affine constraints and provide optimal complexity bounds for four different setups, namely: the function $F(\xb) \triangleq \sum_{i=1}^{m}f_i(\xb)$ is strongly convex and smooth, either strongly convex or smooth or just convex. Our results show that Nesterov's accelerated gradient descent on the dual problem can be executed in a distributed manner and obtains the same optimal rates as in the centralized version of the problem (up to constant or logarithmic factors) with an additional cost related to the spectral gap of the interaction matrix. Finally, we discuss some extensions to the proposed setup such as proximal friendly functions, time-varying graphs, improvement of the condition numbers.