Data-driven linear complexity low-rank approximation of general kernel matrices: A geometric approach

A general, rectangular kernel matrix may be defined as $K_{ij}=\kappa (x_i,y_j)$ where $\kappa (x,y)$ is a kernel function and where $X={\{x_i\}}_{i=1}^m$ and $Y={\{y_i\}}_{i=1}^n$ are two sets of points. In this paper, we seek a low-rank approximation to a kernel matrix where the sets of points $X$ and $Y$ are large and are arbitrarily distributed, such as away from each other, “intermingled”, identical, and so forth. Such rectangular kernel matrices may arise, for example, in Gaussian process regression where $X$ corresponds to the training data and $Y$ corresponds to the test data. In this case, the points are often high-dimensional. Since the point sets are large, we must exploit the fact that the matrix arises from a kernel function, and avoid forming the matrix, and thus ruling out most algebraic techniques. In particular, we seek methods that can scale linearly or nearly linearly with respect to the size of data for a fixed approximation rank. The main idea in this paper is to geometrically select appropriate subsets of points to construct a low rank approximation. An analysis in this paper guides how this selection should be performed.     

Asymptotic spectra of large matrices coming from the symmetrization of Toeplitz structure functions and applications to preconditioning

The singular value distribution of the matrix‐sequence {Y_nT_n[f]}_n , with T_n[f] generated by $f\in L^1([?\pi ,\pi ])$, was shown in [J. Pestana and A.J. Wathen, SIAM J Matrix Anal Appl. 2015;36(1):273‐288]. The results on the spectral distribution of {Y_nT_n[f]}_n were obtained independently in [M. Mazza and J. Pestana, BIT, 59(2):463‐482, 2019] and [P. Ferrari, I. Furci, S. Hon, M.A. Mursaleen, and S. Serra‐Capizzano, SIAM J. Matrix Anal. Appl., 40(3):1066‐1086, 2019]. In the latter reference, the authors prove that {Y_nT_n[f]}_n is distributed in the eigenvalue sense as $${?}_{|f|}(\theta )=\left\{\begin{array}{ll}|f(\theta )|,& \theta \in [0,2\pi ],\\ ?|f(?\theta )|,& \theta \in [?2\pi ,0),\end{array}\right.$$ under the assumptions that f belongs to $L^1([?\pi ,\pi ])$ and has real Fourier coefficients. The purpose of this paper is to extend the latter result to matrix‐sequences of the form {h(T_n[f])}_n , where h is an analytic function. In particular, we provide the singular value distribution of the sequence {h(T_n[f])}_n , the eigenvalue distribution of the sequence {Y_nh(T_n[f])}_n , and the conditions on f and h for these distributions to hold. Finally, the implications of our findings are discussed, in terms of preconditioning and of fast solution methods for the related linear systems.     

A reduced basis approach to large‐scale pseudospectra computation

For studying spectral properties of a nonnormal matrix $A\in {\mathbb{C}}^{n\times n}$, information about its spectrum σ(A) alone is usually not enough. Effects of perturbations on σ(A) can be studied by computing ε‐pseudospectra, i.e. the level sets of the resolvent norm function $\mathsf{g}(z)=\Vert {(zI?A)}^{?1}{\Vert }_2$. The computation of ε‐pseudospectra requires determining the smallest singular values ${\sigma }_{\min }(zI?A)$ for all z on a portion of the complex plane. In this work, we propose a reduced basis approach to pseudospectra computation, which provides highly accurate estimates of pseudospectra in the region of interest, in particular, for pseudospectra estimates in isolated parts of the spectrum containing few eigenvalues of A. It incorporates the sampled singular vectors of zI ? A for different values of z, and implicitly exploits their smoothness properties. It provides rigorous upper and lower bounds for the pseudospectra in the region of interest. In addition, we propose a domain splitting technique for tackling numerically more challenging examples. We present a comparison of our algorithms to several existing approaches on a number of numerical examples, showing that our approach provides significant improvement in terms of computational time.     

Solving the general joint block diagonalization problem via linearly independent eigenvectors of a matrix polynomial

In this paper, we consider the exact/approximate general joint block diagonalization (GJBD) problem of a matrix set ${\{A_i\}}_{i=0}^p$ ( p ≥ 1), where a nonsingular matrix W (often referred to as a diagonalizer) needs to be found such that the matrices W ^HA_iW 's are all exactly/approximately block‐diagonal matrices with as many diagonal blocks as possible. We show that the diagonalizer of the exact GJBD problem can be given by W = [x₁,x₂,…,x_n]Π, where Π is a permutation matrix and x_i's are eigenvectors of the matrix polynomial $P(\lambda )={\sum }_{i=0}^p{\lambda }^i{A}_i$, satisfying that [x₁,x₂,…,x_n] is nonsingular and where the geometric multiplicity of each λ_i corresponding with x_i is equal to 1. In addition, the equivalence of all solutions to the exact GJBD problem is established. Moreover, a theoretical proof is given to show why the approximate GJBD problem can be solved similarly to the exact GJBD problem. Based on the theoretical results, a three‐stage method is proposed, and numerical results show the merits of the method.     

Tensor rank bounds and explicit QTT representations for the inverses of circulant matrices

In this paper, we are concerned with the inversion of circulant matrices and their quantized tensor-train (QTT) structure. In particular, we show that the inverse of a complex circulant matrix $A$, generated by the first column of the form ${(a_0,\dots ,a_{m-1},0,\dots ,0,a_{-n},\dots ,a_{-1})}^{\top }$ admits a QTT representation with the QTT ranks bounded by $(m+n)$. Under certain assumptions on the entries of $A$, we also derive an explicit QTT representation of $A^{-1}$. The latter can be used, for instance, to overcome stability issues arising when numerically solving differential equations with periodic boundary conditions in the QTT format.     

Constructing the field of values of decomposable and general matrices using the ZNN based path following method

This paper describes and develops a fast and accurate path following algorithm that computes the field of values boundary curve $\partial F(A)$ for every conceivable complex or real square matrix $A$. It relies on the matrix flow decomposition algorithm that finds a proper block-diagonal flow representation for the associated hermitean matrix flow ${ℱ}_A(t)=\cos (t)H+\sin (t)K$ under unitary similarity if that is possible. Here ${ℱ}_A(t)$ is the 1-parameter-varying linear combination of the real and skew part matrices $H=(A+A^{\ast })/2$ and $K=(A-A^{\ast })/(2i)$ of $A$. For indecomposable matrix flows, ${ℱ}_A(t)$ has just one block and the ZNN based field of values algorithm works with ${ℱ}_A(t)$ directly. For decomposing flows ${ℱ}_A(t)$, the algorithm decomposes the given matrix $A$ unitarily into block-diagonal form $U^{\ast }AU=\text{diag}(A_j)$ with $j>1$ diagonal blocks $A_j$ whose individual sizes add up to the size of $A$. It then computes the field of values boundaries separately for each diagonal block $A_j$ using the path following ZNN eigenvalue method. The convex hull of all sub-fields of values boundary points $\partial F(A_j)$ finally determines the field of values boundary curve correctly for decomposing matrices $A$. The algorithm removes standard restrictions for path following FoV methods that generally cannot deal with decomposing matrices $A$ due to possible eigencurve crossings of ${ℱ}_A(t)$. Tests and numerical comparisons are included. Our ZNN based method is coded for sequential and parallel computations and both versions run very accurately and fast when compared with Johnson's Francis QR eigenvalue and Bendixson rectangle based method and compute global eigenanalyses of ${ℱ}_A(t_k)$ for large discrete sets of angles $t_k\in [0,2\pi ]$ more slowly.     

Spectral condition‐number estimation of large sparse matrices

We describe a randomized Krylov‐subspace method for estimating the spectral condition number of a real matrix A or indicating that it is numerically rank deficient. The main difficulty in estimating the condition number is the estimation of the smallest singular value ${\sigma }_{\min }$ of A. Our method estimates this value by solving a consistent linear least squares problem with a known solution using a specific Krylov‐subspace method called LSQR. In this method, the forward error tends to concentrate in the direction of a right singular vector corresponding to ${\sigma }_{\min }$. Extensive experiments show that the method is able to estimate well the condition number of a wide array of matrices. It can sometimes estimate the condition number when running dense singular value decomposition would be impractical due to the computational cost or the memory requirements. The method uses very little memory (it inherits this property from LSQR), and it works equally well on square and rectangular matrices.     

Conservative discretizations and parameter‐robust preconditioners for Biot and multiple‐network flux‐based poroelasticity models

The parameters in the governing system of partial differential equations of multiple‐network poroelasticity models typically vary over several orders of magnitude, making its stable discretization and efficient solution a challenging task. In this paper, we prove the uniform Ladyzhenskaya–Babu?ka–Brezzi (LBB) condition and design uniformly stable discretizations and parameter‐robust preconditioners for flux‐based formulations of multiporosity/multipermeability systems. Novel parameter‐matrix‐dependent norms that provide the key for establishing uniform LBB stability of the continuous problem are introduced. As a result, the stability estimates presented here are uniform not only with respect to the Lamé parameter λ but also to all the other model parameters, such as the permeability coefficients K_i; storage coefficients $c_{p_i}$; network transfer coefficients β_{i j},i,j = 1,…,n; the scale of the networks n; and the time step size τ. Moreover, strongly mass‐conservative discretizations that meet the required conditions for parameter‐robust LBB stability are suggested and corresponding optimal error estimates proved. The transfer of the canonical (norm‐equivalent) operator preconditioners from the continuous to the discrete level lays the foundation for optimal and fully robust iterative solution methods. The theoretical results are confirmed in numerical experiments that are motivated by practical applications.     

Large monochromatic components in 3‐edge‐colored Steiner triple systems

It is known that in any r‐coloring of the edges of a complete r‐uniform hypergraph, there exists a spanning monochromatic component. Given a Steiner triple system on n vertices, what is the largest monochromatic component one can guarantee in an arbitrary 3‐coloring of the edges? Gyárfás proved that $\left(2n+3\right)/3$ is an absolute lower bound and that this lower bound is best possible for infinitely many $n$. On the other hand, we prove that for almost all Steiner triple systems the lower bound is actually $\left(1?o\left(1\right)\right)n$. We obtain this result as a consequence of a more general theorem which shows that the lower bound depends on the size of a largest 3‐partite hole (ie, disjoint sets $X_1,X_2,X_3$ with $|X_1|=|X_2|=|X_3|$ such that no edge intersects all of $X_1,X_2,X_3$) in the Steiner triple system (Gyárfás previously observed that the upper bound depends on this parameter). Furthermore, we show that this lower bound is tight unless the structure of the Steiner triple system and the coloring of its edges are restricted in a certain way. We also suggest a variety of other Ramsey problems in the setting of Steiner triple systems.     

Parameter estimation and asymptotic stability in stochastic filtering

In this paper, we study the problem of estimating a Markov chain X$X$ (signal) from its noisy partial information Y$Y$, when the transition probability kernel depends on some unknown parameters. Our goal is to compute the conditional distribution process P{_Xn∣_Yn,…,_Y1}$\mathbb{P}\{X_n\mid Y_n,\dots ,Y_1\}$, referred to hereafter as the optimal filter. Following a standard Bayesian technique, we treat the parameters as a non-dynamic component of the Markov chain. As a result, the new Markov chain is not going to be mixing, even if the original one is. We show that, under certain conditions, the optimal filters are still going to be asymptotically stable with respect to the initial conditions. Thus, by computing the optimal filter of the new system, we can estimate the signal adaptively.     

Asymptotics for the eigenvalues of Toeplitz matrices with a symbol having a power singularity

The present work is devoted to the construction of an asymptotic expansion for the eigenvalues of a Toeplitz matrix $T_n(a)$ as $n$ goes to infinity, with a continuous and real-valued symbol $a$ having a power singularity of degree $\gamma$ with , at one point. The resulting matrix is dense and its entries decrease slowly to zero when moving away from the main diagonal, we apply the so called simple-loop (SL) method for constructing and justifying a uniform asymptotic expansion for all the eigenvalues. Note however, that the considered symbol does not fully satisfy the conditions imposed in previous works, but only in a small neighborhood of the singularity point. In the present work: (i) We construct and justify the asymptotic formulas of the SL method for the eigenvalues ${\lambda }_j(T_n(a))$ with $j⩾\epsilon n$, where the eigenvalues are arranged in nondecreasing order and $\epsilon$ is a sufficiently small fixed number. (ii) We show, with the help of numerical calculations, that the obtained formulas give good approximations in the case . (iii) We numerically show that the main term of the asymptotics for eigenvalues with , formally obtained from the formulas of the SL method, coincides with the main term of the asymptotics constructed and justified in the classical works of Widom and Parter.     

Lower and upper bounds of condition number for Vandermonde-wise matrices and method of fundamental solutions using pseudo radial-lines

Consider the method of fundamental solutions (MFS) for 2D Laplace's equation in a bounded simply connected domain $S$. In the standard MFS, the source nodes are located on a closed contour outside the domain boundary $\mathrm{\Gamma }(=\partial S)$, which is called pseudo-boundary. For circular, elliptic, and general closed pseudo-boundaries, analysis and computation have been studied extensively. New locations of source nodes are proposed along two pseudo radial-lines outside $\mathrm{\Gamma }$. Numerical results are very encouraging and promising. Since the success of the MFS mainly depends on stability, our efforts are focused on deriving the lower and upper bounds of condition number (Cond). The study finds stability properties of new Vandermonde-wise matrices on nodes $x_i\in [a,b]$ with . The Vandermonde-wise matrix is called in this article if it can be decomposed into the standard Vandermonde matrix. New lower and upper bounds of Cond are first derived for the standard Vandermonde matrix, and then for new algorithms of the MFS using two pseudo radial-lines. Both lower and upper bounds of Cond are intriguing in the stability study for the MFS. Numerical experiments are carried out to verify the stability analysis made. Since the fundamental solutions (as $\{\ln |\overline{{PQ}_i}|\}$) are the basis functions of the MFS, new Vandermonde-wise matrices are found. Since the nodes $x_i\in [a,b]$ with may come from approximations and interpolations by the Laurent polynomials with singular part, the conclusions in this article are important not only to the MFS but also to matrix analysis.     

Spectral analysis of Pk Finite Element matrices in the case of Friedrichs–Keller triangulations via Generalized Locally Toeplitz technology

In the present article, we consider a class of elliptic partial differential equations with Dirichlet boundary conditions and where the operator is div(?a( x )?·), with a continuous and positive over $\stackrel{￣}{\mathrm{\Omega }}$, Ω being an open and bounded subset of ${\mathbb{R}}^d$, d≥1. For the numerical approximation, we consider the classical ${\mathbb{P}}_k$ Finite Elements, in the case of Friedrichs–Keller triangulations, leading, as usual, to sequences of matrices of increasing size. The new results concern the spectral analysis of the resulting matrix‐sequences in the direction of the global distribution in the Weyl sense, with a concise overview on localization, clustering, extremal eigenvalues, and asymptotic conditioning. We study in detail the case of constant coefficients on Ω=(0,1)² and we give a brief account in the more involved case of variable coefficients and more general domains. Tools are drawn from the Toeplitz technology and from the rather new theory of Generalized Locally Toeplitz sequences. Numerical results are shown for a practical evidence of the theoretical findings.     

On the fixed part of pluricanonical systems for surfaces

We show that $|m{K}_X|$ defines a birational map and has no fixed part for some bounded positive integer m for any $\frac{1}{2}$-lc surface X such that $K_X$ is big and nef. For every positive integer $n\ge 3$, we construct a sequence of projective surfaces $X_{n,i}$, such that $K_{X_{n,i}}$ is ample, $\mathrm{mld}(X_{n,i})>\frac{1}{n}$ for every i, ${\lim }_{i\to +\infty }\mathrm{mld}(X_{n,i})=\frac{1}{n}$, and for any positive integer m, there exists i such that $|m{K}_{X_{n,i}}|$ has nonzero fixed part. These results answer the surface case of a question of Xu.     

Maximal ℓp‐regularity of multiterm fractional equations with delay

We provide a characterization for the existence and uniqueness of solutions in the space of vector‐valued sequences ${?}_p(?,X)$ for the multiterm fractional delayed model in the form $${\mathrm{\Delta }}^{\alpha }u(n)+\lambda {\mathrm{\Delta }}^{\beta }u(n)=Au(n)+u(n?\tau )+f(n),n\in ?,\alpha ,\beta \in {?}_{+},\tau \in ?,\lambda \in ?,$$ where X is a Banach space, A is a closed linear operator with domain D(A) defined on X, $f\in {?}_p(?,X)$ and Δ^Γ denotes the Grünwald–Letkinov fractional derivative of order Γ > 0. We also give some conditions to ensure the existence of solutions when adding nonlinearities. Finally, we illustrate our results with an example given by a general abstract nonlinear model that includes the fractional Fisher equation with delay.     

Fast solution of three-dimensional elliptic equations with randomly generated jumping coefficients by using tensor-structured preconditioners

In this paper, we propose and analyze the numerical algorithms for fast solution of periodic elliptic problems in random media in ${ℝ}^d$, $d=2,3$. Both the two-dimensional (2D) and three-dimensional (3D) elliptic problems are considered for the jumping equation coefficients built as a checkerboard type configuration of bumps randomly distributed on a large $L\times L$, or $L\times L\times L$ lattice, respectively. The finite element method discretization procedure on a 3D $n\times n\times n$ uniform tensor grid is described in detail, and the Kronecker tensor product approach is proposed for fast generation of the stiffness matrix. We introduce tensor techniques for the construction of the low Kronecker rank spectrally equivalent preconditioner in a periodic setting to be used in the framework of the preconditioned conjugate gradient iteration. The discrete 3D periodic Laplacian pseudo-inverse is first diagonalized in the Fourier basis, and then the diagonal matrix is reshaped into a fully populated third-order tensor of size $n\times n\times n$. The latter is approximated by a low-rank canonical tensor by using the multigrid Tucker-to-canonical tensor transform. As an example, we apply the presented solver in numerical analysis of stochastic homogenization method where the 3D elliptic equation should be solved many hundred times, and where for every random sampling of the equation coefficient one has to construct the new stiffness matrix and the right-hand side. The computational characteristics of the presented solver in terms of a lattice parameter $L$ and the grid-size, $n^d$, in both 2D and 3D cases are illustrated in numerical tests. Our solver can be used in various applications where the elliptic problem should be solved for a number of different coefficients for example, in many-particle dynamics, protein docking problems or stochastic modeling.     

Strong Asymptotics of Planar Orthogonal Polynomials: Gaussian Weight Perturbed by Finite Number of Point Charges

We consider the orthogonal polynomial pn(z) with respect to the planar measure supported on the whole complex plane $${\mathrm{e}}^{-N|z{|}^2}\prod _{j=1}^{\nu }|z-a_j{|}^{2c_j}\mathrm{d}A(z)$$ where dA is the Lebesgue measure of the plane, N is a positive constant, {c₁, …, c_ν} are nonzero real numbers greater than −1 and $\{a_1,\text{…},a_{\nu }\}\subset \mathbb{D}\setminus \{0\}$ are distinct points inside the unit disk. In the scaling limit when n/N = 1 and n → ∞ we obtain the strong asymptotics of the polynomial p_n(z). We show that the support of the roots converges to what we call the “multiple Szegő curve,” a certain connected curve having ν + 1 components in its complement. We apply the nonlinear steepest descent method [9,10] on the matrix Riemann-Hilbert problem of size (ν + 1) × (ν + 1) posed in [22]. © 2023 Wiley Periodicals, LLC.     

Theorems of Wolff–Denjoy type for semigroups of nonexpansive mappings in geodesic spaces

We show that if $S={\{f_t:Y\to Y\}}_{t\ge 0}$ is a one-parameter continuous semigroup of nonexpansive mappings acting on a complete locally compact geodesic space $(Y,d)$ that satisfies some geometric properties, then there exists $\xi \in \partial Y$ such that S converge uniformly on bounded sets of Y to ξ. In particular, our result applies to strictly convex bounded domains in ${\mathbb{R}}^n$ or ${\mathbb{C}}^n$ with respect to a large class of metrics including Hilbert's and Kobayashi's metrics.     

Semi‐active ℋ∞ damping optimization by adaptive interpolation

In this work we consider the problem of semi‐active damping optimization of mechanical systems with fixed damper positions. Our goal is to compute a damping that is locally optimal with respect to the ${?}_{\infty }$‐norm of the transfer function from the exogenous inputs to the performance outputs. We make use of a new greedy method for computing the ${?}_{\infty }$‐norm of a transfer function based on rational interpolation. In this paper, this approach is adapted to parameter‐dependent transfer functions. The interpolation leads to parametric reduced‐order models that can be optimized more efficiently. At the optimizers we then take new interpolation points to refine the reduced‐order model and to obtain updated optimizers. In our numerical examples we show that this approach normally converges fast and thus can highly accelerate the optimization procedure. Another contribution of this work is heuristics for choosing initial interpolation points.     

I‐convexity and Q‐convexity in Orlicz–Bochner function spaces endowed with the Orlicz norm

Some characterizations of I‐convexity and Q‐convexity of Banach space are obtained. Moreover, the criteria is shown for Orlicz–Bochner function spaces $L_M(\mu ,X)$ endowed with the Orlicz norm being I‐convex as well as being Q‐convex.