Search

Let $A$ be a square matrix with a given structure (e.g. real matrix, sparsity pattern, Toeplitz structure, etc.) and assume that it is unstable, i.e. at least one of its eigenvalues lies in the complex right half-plane. The problem of stabilizing $A$ consists in the computation of a matrix $B$, whose eigenvalues have negative real part and such that the perturbation $\Delta=B-A$ has minimal norm. The structured stabilization further requires that the perturbation preserves the structural pattern of $A$. We solve this non-convex problem by a two-level procedure which involves the computation of the stationary points of a matrix ODE. We exploit the low rank underlying features of the problem by using an adaptive-rank integrator that follows slavishly the rank of the solution. We show the benefits derived from the low rank setting in several numerical examples, which also allow to deal with high dimensional problems., Comment: 17 pages, 5 figures, 7 tables

We propose a novel methodology to solve a key eigenvalue optimization problem which arises in the contractivity analysis of neural ODEs. When looking at contractivity properties of a one layer weight-tied neural ODE $\dot{u}(t)=\sigma(Au(t)+b)$ (with $u,b \in {\mathbb R}^n$, $A$ is a given $n \times n$ matrix, $\sigma : {\mathbb R} \to {\mathbb R}^+$ denotes an activation function and for a vector $z \in {\mathbb R}^n$, $\sigma(z) \in {\mathbb R}^n$ has to be interpreted entry-wise), we are led to study the logarithmic norm of a set of products of type $D A$, where $D$ is a diagonal matrix such that ${\mathrm{diag}}(D) \in \sigma'({\mathbb R}^n)$. Specifically, given a real number $c$ (usually $c=0$), the problem consists in finding the largest positive interval $\chi\subseteq \mathbb [0,\infty)$ such that the logarithmic norm $\mu(DA) \le c$ for all diagonal matrices $D$ with $D_{ii}\in \chi$. We propose a two-level nested methodology: an inner level where, for a given $\chi$, we compute an optimizer $D^\star(\chi)$ by a gradient system approach, and an outer level where we tune $\chi$ so that the value $c$ is reached by $\mu(D^\star(\chi)A)$. We extend the proposed two-level approach to the general multilayer, and possibly time-dependent, case $\dot{u}(t) = \sigma( A_k(t) \ldots \sigma ( A_{1}(t) u(t) + b_{1}(t) ) \ldots + b_{k}(t) )$ and we propose several numerical examples to illustrate its behaviour, including its stabilizing performance on a one-layer neural ODE applied to the classification of the MNIST handwritten digits dataset., Comment: 23 pages, 5 figures, 3 tables

Recent work in deep learning has shown strong empirical and theoretical evidence of an implicit low-rank bias: weight matrices in deep networks tend to be approximately low-rank and removing relatively small singular values during training or from available trained models may significantly reduce model size while maintaining or even improving model performance. However, the majority of the theoretical investigations around low-rank bias in neural networks deal with oversimplified deep linear networks. In this work, we consider general networks with nonlinear activations and the weight decay parameter, and we show the presence of an intriguing neural rank collapse phenomenon, connecting the low-rank bias of trained networks with networks' neural collapse properties: as the weight decay parameter grows, the rank of each layer in the network decreases proportionally to the within-class variability of the hidden-space embeddings of the previous layers. Our theoretical findings are supported by a range of experimental evaluations illustrating the phenomenon.

Techniques based on $k$-th order Hodge Laplacian operators $L_k$ are widely used to describe the topology as well as the governing dynamics of high-order systems modeled as simplicial complexes. In all of them, it is required to solve a number of least square problems with $L_k$ as coefficient matrix, for example in order to compute some portions of the spectrum or integrate the dynamical system. In this work, we introduce the notion of optimal collapsible subcomplex and we present a fast combinatorial algorithm for the computation of a sparse Cholesky-like preconditioner for $L_k$ that exploits the topological structure of the simplicial complex. The performance of the preconditioner is tested for conjugate gradient method for least square problems (CGLS) on a variety of simplicial complexes with different dimensions and edge densities. We show that, for sparse simplicial complexes, the new preconditioner reduces significantly the condition number of $L_k$ and performs better than the standard incomplete Cholesky factorization., Comment: 22 pages, 8 figures

There exist excellent codes for an efficient numerical treatment of stiff and differential-algebraic problems. Let us mention {\sc Radau5} which is based on the $3$-stage Radau IIA collocation method, and its extension to problems with discrete delays {\sc Radar5}. The aim of the present work is to present a technique that permits a direct application of these codes to problems having a right-hand side with an additional distributed delay term (which is a special case of an integro-differential equation). Models with distributed delays are of increasing importance in pharmacodynamics and pharmacokinetics for the study of the interaction between drugs and the body. The main idea is to approximate the distribution kernel of the integral term by a sum of exponential functions or by a quasi-polynomial expansion, and then to transform the distributed (integral) delay term into a set of ordinary differential equations. This set is typically stiff and, for some distribution kernels (e.g., Pareto distribution), it contains discrete delay terms with constant delay. The original equations augmented by this set of ordinary differential equations can have a very large dimension, and a careful treatment of the solution of the arising linear systems is necessary. The use of the codes {\sc Radau5} and {\sc Radar5} is illustrated at three examples (two test equations and one problem taken from pharmacodynamics). The driver programs for these examples are publicly available from the homepages of the authors., Comment: 20 pages, 3 figures

Given a matrix-valued function $\mathcal{F}(\lambda)=\sum_{i=1}^d f_i(\lambda) A_i$, with complex matrices $A_i$ and $f_i(\lambda)$ entire functions for $i=1,\ldots,d$, we discuss a method for the numerical approximation of the distance to singularity of $\mathcal{F}(\lambda)$. The closest singular matrix-valued function $\widetilde{\mathcal{F}}(\lambda)$ with respect to the Frobenius norm is approximated using an iterative method. The property of singularity on the matrix-valued function is translated into a numerical constraint for a suitable minimization problem. Unlike the case of matrix polynomials, in the general setting of matrix-valued functions the main issue is that the function $\det ( \widetilde{\mathcal{F}}(\lambda) )$ may have an infinite number of roots. An important feature of the numerical method consists in the possibility of addressing different structures, such as sparsity patterns induced by the matrix coefficients, in which case the search of the closest singular function is restricted to the class of functions preserving the structure of the matrices., Comment: 24 pages, 5 figures

In a recent paper by Guglielmi and Hairer (SIADS 2015), an analysis in the $\varepsilon\to 0$ limit was proposed of regularized discontinuous ODEs in codimension-2 switching domains; this was obtained by studying a certain 2-dimensional system describing the so-called hidden dynamics. In particular, the existence of a unique limit solution was not proved in all cases, a few of which were labeled as ambiguous, and it was not clear whether or not the ambiguity could be resolved. In this paper, we show that it cannot be resolved in general. A first contribution of this paper is an illustration of the dependence of the limit solution on the form of the switching function. Considering the parameter dependence in the ambiguous class of discontinuous systems, a second contribution is a bifurcation analysis, revealing a range of possible behaviors. Finally, we investigate the sensitivity of solutions in the transition from codimension-2 domains to codimension-3 when there is a limit cycle in the hidden dynamics., Comment: 32 pages, 5 figures

Spectral clustering is a well-known technique which identifies $k$ clusters in an undirected graph with weight matrix $W\in\mathbb{R}^{n\times n}$ by exploiting its graph Laplacian $L(W)$, whose eigenvalues $0=\lambda_1\leq \lambda_2 \leq \dots \leq \lambda_n$ and eigenvectors are related to the $k$ clusters. Since the computation of $\lambda_{k+1}$ and $\lambda_k$ affects the reliability of this method, the $k$-th spectral gap $\lambda_{k+1}-\lambda_k$ is often considered as a stability indicator. This difference can be seen as an unstructured distance between $L(W)$ and an arbitrary symmetric matrix $L_\star$ with vanishing $k$-th spectral gap. A more appropriate structured distance to ambiguity such that $L_\star$ represents the Laplacian of a graph has been proposed by Andreotti et al. (2021). Slightly differently, we consider the objective functional $ F(\Delta)=\lambda_{k+1}\left(L(W+\Delta)\right)-\lambda_k\left(L(W+\Delta)\right)$, where $\Delta$ is a perturbation such that $W+\Delta$ has non-negative entries and the same pattern of $W$. We look for an admissible perturbation $\Delta_\star$ of smallest Frobenius norm such that $F(\Delta_\star)=0$. In order to solve this optimization problem, we exploit its low rank underlying structure. We formulate a rank-4 symmetric matrix ODE whose stationary points are the optimizers sought. The integration of this equation benefits from the low rank structure with a moderate computational effort and memory requirement, as it is shown in some illustrative numerical examples., Comment: 24 pages, 4 figures, 5 tables

Consider a matrix polynomial $P \left( \lambda \right)= A_0 + \lambda A_1 + \ldots + \lambda^d A_d$, with $A_0,\ldots, A_d$ complex (or real) matrices with a certain structure. In this paper we discuss an iterative method to numerically approximate the closest structured singular matrix polynomial $\widetilde P\left( \lambda \right)$, using the distance induced by the Frobenius norm. An important peculiarity of the approach we propose is the possibility to include different types of structural constraints. The method also allows us to limit the perturbations to just a few matrices and also to include additional structures, such as the preservation of the sparsity pattern of one or more matrices $A_i$, and also collective-like properties, like a palindromic structure. The iterative method is based on the numerical integration of the gradient system associated with a suitable functional which quantifies the distance to singularity of a matrix polynomial., Comment: 28 pages

The homology groups of a simplicial complex reveal fundamental properties of the topology of the data or the system and the notion of topological stability naturally poses an important yet not fully investigated question. In the current work, we study the stability in terms of the smallest perturbation sufficient to change the dimensionality of the corresponding homology group. Such definition requires an appropriate weighting and normalizing procedure for the boundary operators acting on the Hodge algebra's homology groups. Using the resulting boundary operators, we then formulate the question of structural stability as a spectral matrix nearness problem for the corresponding higher-order graph Laplacian. We develop a bilevel optimization procedure suitable for the formulated matrix nearness problem and illustrate the method's performance on a variety of synthetic quasi-triangulation datasets and transportation networks., Comment: 25 pages, 9 figures

A new approach to solving eigenvalue optimization problems for large structured matrices is proposed and studied. The class of optimization problems considered is related to computing structured pseudospectra and their extremal points, and to structured matrix nearness problems such as computing the structured distance to instability or to singularity. The structure can be a general linear structure and includes, for example, large matrices with a given sparsity pattern, matrices with given range and co-range, and Hamiltonian matrices. Remarkably, the eigenvalue optimization can be performed on the manifold of complex (or real) rank-1 matrices, which yields a significant reduction of storage and in some cases of the computational cost. The method relies on a constrained gradient system and the projection of the gradient onto the tangent space of the manifold of complex rank-$1$ matrices. It is shown that near a local minimizer this projection is very close to the identity map, and so the computationally favorable rank-1 projected system behaves locally like the %computationally expensive gradient system.

In this paper we discuss a projection model order reduction (MOR) method for a class of parametric linear evolution PDEs, which is based on the application of the Laplace transform. The main advantage of this approach consists in the fact that, differently from time stepping methods, like Runge-Kutta integrators, the Laplace transform allows to compute the solution directly at a given instant, which can be done by approximating the contour integral associated to the inverse Laplace transform by a suitable quadrature formula. In terms of some classical MOR methodology, this determines a significant improvement in the reduction phase - like the one based on the classical proper orthogonal decomposition (POD) - since the number of vectors to which the decomposition applies is drastically reduced as it does not contain all intermediate solutions generated along an integration grid by a time stepping method. We show the effectiveness of the method by some illustrative parabolic PDEs arising from finance and also provide some evidence that the method we propose, when applied to a linear advection equation, does not suffer the problem of slow decay of singular values which instead affects time stepping methods for the numerical approximation of the Cauchy problem arising from space discretization.

COVID-19 data released by public health authorities features the presence of notable time-delays, corresponding to the difference between actual time of infection and identification of infection. These delays have several causes, including the natural incubation period of the virus, availability and speed of testing facilities, population demographics, and testing center capacity, among others. Such delays have important ramifications for both the mathematical modeling of COVID-19 contagion and the design and evaluation of intervention strategies. In the present work, we introduce a novel optimization technique for the identification of time delays in COVID-19 data, making use of a delay-differential equation model. The proposed method is general in nature and may be applied not only to COVID-19, but for generic dynamical systems in which time delays may be present.

We study computational methods for computing the distance to singularity, the distance to the nearest high index problem, and the distance to instability for linear differential-algebraic systems (DAEs) with dissipative Hamiltonian structure. While for general unstructured DAEs the characterization of these distances is very difficult, and partially open, it has been recently shown that for dissipative Hamiltonian systems and related matrix pencils there exist explicit characterizations. We will use these characterizations for the development of computational methods to compute these distances via methods that follow the flow of a differential equation converging to the smallest perturbation that destroys the property of regularity, index one or stability.

In the wake of the 2020 COVID-19 epidemic, much work has been performed on the development of mathematical models for the simulation of the epidemic, and of disease models generally. Most works follow the susceptible-infected-removed (SIR) compartmental framework, modeling the epidemic with a system of ordinary differential equations. Alternative formulations using a partial differential equation (PDE) to incorporate both spatial and temporal resolution have also been introduced, with their numerical results showing potentially powerful descriptive and predictive capacity. In the present work, we introduce a new variation to such models by using delay differential equations (DDEs). The dynamics of many infectious diseases, including COVID-19, exhibit delays due to incubation periods and related phenomena. Accordingly, DDE models allow for a natural representation of the problem dynamics, in addition to offering advantages in terms of computational time and modeling, as they eliminate the need for additional, difficult-to-estimate, compartments (such as exposed individuals) to incorporate time delays. Here, we introduce a DDE epidemic model in both an ordinary- and partial differential equation framework. We present a series of mathematical results assessing the stability of the formulation. We then perform several numerical experiments, validating both the mathematical results and establishing model's ability to reproduce measured data on realistic problems.

We generalize ideas in the recent literature and develop new ones in order to propose a general class of contour integral methods for linear convection-diffusion PDEs and in particular for those arising in finance. These methods aim to provide a numerical approximation of the solution by computing its inverse Laplace transform. The choice of the integration contour is determined by the computation of a few suitably weighted pseudo-spectral level sets of the leading operator of the equation. Parabolic and hyperbolic profiles proposed in the literature are investigated and compared to the elliptic contour originally proposed by Guglielmi, L\'opez-Fern\'andez and Nino. In summary, the article (i) provides a comparison among three different integration profiles; (ii) proposes a new fast pseudospectral roaming method; (iii) optimizes the selection of time windows on which one may arbitrarily approximate the solution by no extra computational cost with respect to the case of a fixed time instant; (iv) focuses extensively on computational aspects and it is the reference of the MATLAB code https://github.com/MattiaManucci/Contour_Integral_Methods.git, where all algorithms described here are implemented.

We propose and study an algorithm for computing a nearest passive system to a given non-passive linear time-invariant system (with much freedom in the choice of the metric defining `nearest', which may be restricted to structured perturbations), and also a closely related algorithm for computing the structured distance of a given passive system to non-passivity. Both problems are addressed by solving eigenvalue optimization problems for Hamiltonian matrices that are constructed from perturbed system matrices. The proposed algorithms are two-level methods that optimize the Hamiltonian eigenvalue of smallest positive real part over perturbations of a fixed size in the inner iteration, using a constrained gradient flow. They optimize over the perturbation size in the outer iteration, which is shown to converge quadratically in the typical case of a defective coalescence of simple eigenvalues approaching the imaginary axis. For large systems, we propose a variant of the algorithm that takes advantage of the inherent low-rank structure of the problem. Numerical experiments illustrate the behavior of the proposed algorithms.

Rank deficient Hankel matrices are at the core of several applications. However, in practice, the coefficients of these matrices are noisy due to e.g. measurements errors and computational errors, so generically the involved matrices are full rank. This motivates the problem of Hankel structured low-rank approximation. Structured low-rank approximation problems, in general, do not have a global and efficient solution technique. In this paper we propose a local optimization approach based on a two-levels iteration. Experimental results show that the proposed algorithm usually achieves good accuracy and shows a higher robustness with respect to the initial approximation, compared to alternative approaches.

We study asymptotic stability of continuous-time systems with mode-dependent guaranteed dwell time. These systems are reformulated as special cases of a general class of mixed (discrete-continuous) linear switching systems on graphs, in which some modes correspond to discrete actions and some others correspond to continuous-time evolutions. Each discrete action has its own positive weight which accounts for its time-duration. We develop a theory of stability for the mixed systems; in particular, we prove the existence of an invariant Lyapunov norm for mixed systems on graphs and study its structure in various cases, including discrete-time systems for which discrete actions have inhomogeneous time durations. This allows us to adapt recent methods for the joint spectral radius computation (Gripenberg's algorithm and the Invariant Polytope Algorithm) to compute the Lyapunov exponent of mixed systems on graphs.

As an indicator of the stability of spectral clustering of an undirected weighted graph into $k$ clusters, the $k$th spectral gap of the graph Laplacian is often considered. The $k$th spectral gap is characterized in this paper as an unstructured distance to ambiguity, namely as the minimal distance of the Laplacian to arbitrary symmetric matrices with vanishing $k$th spectral gap. As a conceptually more appropriate measure of stability, the structured distance to ambiguity of the $k$-clustering is introduced as the minimal distance of the Laplacian to Laplacians of graphs with the same vertices and edges but with weights that are perturbed such that the $k$th spectral gap vanishes. To compute a solution to this matrix nearness problem, a two-level iterative algorithm is proposed that uses a constrained gradient system of matrix differential equations in the inner iteration and a one-dimensional optimization of the perturbation size in the outer iteration. The structured and unstructured distances to ambiguity are compared on some example graphs. The numerical experiments show, in particular, that selecting the number $k$ of clusters according to the criterion of maximal stability can lead to different results for the structured and unstructured stability indicators.

Given a matrix $A\in \R^{n\times n}$ and a tall rectangular matrix $B \in \R^{n\times q}$, $q < n$, we consider the problem of making the pair $(A,B)$ dissipative, that is the determination of a {\it feedback} matrix $K \in \R^{q\times n}$ such that the field of values of $A-B K$ lies in the left half open complex plane. We review and expand classical results available in the literature on the existence and parameterization of the class of dissipating matrices, and we explore new matrix properties associated with the problem. In addition, we discuss various computational strategies for approximating the minimal Frobenius norm dissipating $K$.

We consider the problem of computing the closest stable/unstable non-negative matrix to a given real matrix. This problem is important in the study of linear dynamical systems, numerical methods, etc. The distance between matrices is measured in the Frobenius norm. The problem is addressed for two types of stability: the Schur stability (the matrix is stable if its spectral radius is smaller than one) and the Hurwitz stability (the matrix is stable if its spectral abscissa is negative). We show that the closest unstable matrix can always be explicitly found. For the closest stable matrix, we present an iterative algorithm which converges to a local minimum with a linear rate. It is shown that the total number of local minima can be exponential in the dimension. Numerical results and the complexity estimates are presented.

In this paper a novel contour integral method is proposed for linear convection-diffusion equations. The method is based on the inversion of the Laplace transform and makes use of a contour given by an elliptic arc joined symmetrically to two half-lines. The trapezoidal rule is the chosen integration method for the numerical inversion of the Laplace transform, due to its well-known fast convergence properties when applied to analytic functions. Error estimates are provided as well as careful indications about the choice of several involved parameters. The method selects the elliptic arc in the integration contour by an algorithmic strategy based on the computation of pseudospectral level sets of the discretized differential operator. In this sense the method is general and can be applied to any linear convection-diffusion equation without knowing any a priori information about its pseudospectral geometry. Numerical experiments performed on the Black-Scholes ($1D$) and Heston ($2D$) equations show that the method is competitive with other contour integral methods available in the literature.

Given a connected undirected weighted graph, we are concerned with problems related to partitioning the graph. First of all we look for the closest disconnected graph (the minimum cut problem), here with respect to the Euclidean norm. We are interested in the case of constrained minimum cut problems, where constraints include cardinality or membership requirements, which leads to NP-hard combinatorial optimization problems. Furthermore, we are interested in ambiguity issues, that is in the robustness of clustering algorithms that are based on Fiedler spectral partitioning. The above-mentioned problems are restated as matrix nearness problems for the weight matrix of the graph. A key element in the solution of these matrix nearness problems is the use of a constrained gradient system of matrix differential equations.

We present several results describing the interplay between the max algebraic joint spectral radius (JSR) for compact sets of matrices and suitably defined matrix norms. In particular, we extend a classical result for the conventional algebra, showing that the JSR can be described in terms of induced norms of the matrices in the set. We also show that for a set generating an irreducible semigroup (in a cone-theoretic sense), a monotone Barabanov norm always exists. This fact is then used to show that the max algebraic JSR is locally Lipschitz continuous on the space of compact irreducible sets of matrices with respect to the Hausdorff distance. We then prove that the JSR is Hoelder continuous on the space of compact sets of nonnegative matrices. Finally, we prove a strict monotonicity property for the max algebraic JSR that echoes a fact for the classical JSR.

We propose a fast method to approximate the real stability radius of a linear dynamical system with output feedback, where the perturbations are restricted to be real valued and bounded with respect to the Frobenius norm. Our work builds on a number of scalable algorithms that have been proposed in recent years, ranging from methods that approximate the complex or real pseudospectral abscissa and radius of large sparse matrices (and generalizations of these methods for pseudospectra to spectral value sets) to algorithms for approximating the complex stability radius (the reciprocal of the $H_\infty$ norm). Although our algorithm is guaranteed to find only upper bounds to the real stability radius, it seems quite effective in practice. As far as we know, this is the first algorithm that addresses the Frobenius-norm version of this problem. Because the cost mainly consists of computing the eigenvalue with maximal real part for continuous-time systems (or modulus for discrete-time systems) of a sequence of matrices, our algorithm remains very efficient for large-scale systems provided that the system matrices are sparse.

We consider linear dynamical systems with a structure of a multigraph. The vertices are associated to linear spaces and the edges correspond to linear maps between those spaces. We analyse the asymptotic growth of trajectories (associated to paths along the multigraph), the stability and the stabilizability problems. This generalizes the classical linear switching systems and their recent extensions to Markovian systems, to systems generated by regular languages, etc. We show that an arbitrary system can be factorized into several irreducible systems on strongly connected multigraphs. For the latter systems, we prove the existence of invariant (Barabanov) multinorm and derive a method of its construction. The method works for a vast majority of systems and finds the joint spectral radius (Lyapunov exponent). Numerical examples are presented and applications to the study of fractals, attractors, and multistep methods for ODEs are discussed., Comment: 45 pages

A novel method for approximating structured singular values (also known as mu-values) is proposed and investigated. These quantities constitute an important tool in the stability analysis of uncertain linear control systems as well as in structured eigenvalue perturbation theory. Our approach consists of an inner-outer iteration. In the outer iteration, a Newton method is used to adjust the perturbation level. The inner iteration solves a gradient system associated with an optimization problem on the manifold induced by the structure. Numerical results and comparison with the well-known Matlab function mussv, implemented in the Matlab Control Toolbox, illustrate the behavior of the method.

We consider the optimal control design problem for discrete-time LTI systems with state feedback, when the actuation signal is subject to unmeasurable switching propagation delays, due to e.g. the routing in a multi-hop communication network and/or jitter. In particular, we set up a constrained optimization problem where the cost function is the worst-case $\mathcal{L}_2$ norm for all admissible switching delays. We first show how to model these systems as pure switching linear systems, and as main contribution of the paper we provide an algorithm to compute a sub-optimal solution., Comment: arXiv admin note: text overlap with arXiv:1401.1673

In this paper, we present a new matrix approach for the analysis of subdivision schemes whose non-stationarity is due to linear dependency on parameters whose values vary in a compact set. Indeed, we show how to check the convergence in $C^{\ell}(\RR^s)$ and determine the H\"older regularity of such level and parameter dependent schemes efficiently via the joint spectral radius approach. The efficiency of this method and the important role of the parameter dependency are demonstrated on several examples of subdivision schemes whose properties improve the properties of the corresponding stationary schemes. Moreover, we derive necessary criteria for a function to be generated by some level dependent scheme and, thus, expose the limitations of such schemes.

We generalize the recent invariant polytope algorithm for computing the joint spectral radius and extend it to a wider class of matrix sets. This, in particular, makes the algorithm applicable to sets of matrices that have finitely many spectrum maximizing products. A criterion of convergence of the algorithm is proved. As an application we solve two challenging computational open problems. First we find the regularity of the Butterfly subdivision scheme for various parameters $\omega$. In the "most regular" case $\omega = \frac{1}{16}$, we prove that the limit function has H\"older exponent $2$ and its derivative is "almost Lipschitz" with logarithmic factor $2$. Second we compute the H\"older exponent of Daubechies wavelets of high order., Comment: 36 pages

In this paper, we study scalar multivariate non-stationary subdivision schemes with integer dilation matrix M=mI, m >=2, and present a general approach for checking their convergence and for determining their H\"older regularity. The combination of the concepts of asymptotic similarity and approximate sum rules allows us to link stationary and non-stationary settings and to employ recent advances in methods for exact computation of the joint spectral radius. As an application, we prove a recent conjecture on the H\"older regularity of the generalized Daubechies wavelets. We illustrate our results with several examples.

We present a new approach for constructing polytope Lyapunov functions for continuous-time linear switching systems (LSS). This allows us to decide the stability of LSS and to compute the Lyapunov exponent with a good precision in relatively high dimensions. The same technique is also extended for stabilizability of positive systems by evaluating a polytope concave Lyapunov function ("antinorm") in the cone. The method is based on a suitable discretization of the underlying continuous system and provides both a lower and an upper bound for the Lyapunov exponent. The absolute error in the Lyapunov exponent computation is estimated from above and proved to be linear in the dwell time. The practical efficiency of the new method is demonstrated in several examples and in the list of numerical experiments with randomly generated matrices of dimensions up to $10$ (for general linear systems) and up to $100$ (for positive systems). The development of the method is based on several theoretical results proved in the paper: the existence of monotone invariant norms and antinorms for positively irreducible systems, the equivalence of all contractive norms for stable systems and the linear convergence theorem.

Let $A$ be either a complex or real matrix with all distinct eigenvalues. We propose a new method for the computation of both the unstructured and the real-structured (if the matrix is real) distance $w_{\mathbb K}(A)$ (where ${\mathbb K}=\mathbb C$ if general complex matrices are considered and ${\mathbb K} ={\mathbb R}$ if only real matrices are allowed) of the matrix $A$ from the set of defective matrices, that is the set of those matrices with at least a multiple eigenvalue with algebraic multiplicity larger than its geometric multiplicity. For $0 < \varepsilon \le w_{\mathbb K}(A)$, this problem is closely related to the computation of the most ill-conditioned $\varepsilon$-pseudoeigenvalues of $A$, that is points in the $\varepsilon$-pseudospectrum of $A$ characterized by the highest condition number. The method we propose couples a system of differential equations on a low rank (possibly structured) manifold which computes the $\varepsilon$-pseudoeigenvalue of $A$ which is closest to coalesce, with a fast Newton-like iteration aiming to determine the minimal value $\varepsilon$ such that such an $\varepsilon$-pseudoeigenvalue becomes defective. The method has a local behaviour; this means that in general we find upper bounds for $w_{\mathbb K}(A)$. However, they usually provide good approximations, in those (simple) cases where we can check this. The methodology can be extended to a structured matrix where it is required that the distance is computed within some manifold defining the structure of the matrix. In this paper we extensively examine the case of real matrices but we also consider pattern structures. As far as we know there do not exist methods in the literature able to compute such distance., Comment: 30 pages

We describe new methods for deciding the stability of switching systems. The methods build on two ideas previously appeared in the literature: the polytope norm iterative construction, and the lifting procedure. Moreover, the combination of these two ideas allows us to introduce a pruning algorithm which can importantly reduce the computational burden. We prove several appealing theoretical properties of our methods like a finiteness computational result which extends a known result for unlifted sets of matrices, and provide numerical examples of their good behaviour.