On the numerical approximation of the distance to singularity for matrix-valued functions

TL;DR

The paper addresses the challenge of approximating the distance to singularity for general matrix-valued functions $\mathcal{F}(\lambda)=\sum_{i=1}^d f_i(\lambda) A_i$, where $\det(\mathcal{F}(\lambda))$ may have infinitely many zeros. It introduces a two-level, ODE-based optimization method that discretizes the problem on the unit circle via $m$ sampling points $\mu_j$, turning the continuous singularity constraint into a finite set of conditions enforced through the smallest singular values at those points; an inner gradient-flow minimizes a sum of squared singular values, while an outer Newton–bisection step recovers the minimal perturbation size. The framework accommodates structured perturbations through projections onto subspaces, and extends naturally to matrix polynomials with improved point-count efficiency, yielding practical benefits for stability analysis in delay systems and nonlinear eigenvalue problems. The numerical experiments demonstrate effectiveness, robustness, and speed-ups in both unstructured and structured settings, including comparisons to existing polynomial-specific methods and applicability to delay-differential equation contexts.

Abstract

Given a matrix-valued function $\mathcal{F}(λ)=\sum_{i=1}^d f_i(λ) A_i$, with complex matrices $A_i$ and $f_i(λ)$ entire functions for $i=1,\ldots,d$, we discuss a method for the numerical approximation of the distance to singularity of $\mathcal{F}(λ)$. The closest singular matrix-valued function $\widetilde{\mathcal{F}}(λ)$ 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}}(λ) )$ 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.

Figures (7)

• Figure 1: Plot of the solution of the original problem (left) and error wrt perturbed problem (right) in the case of delay $\tau=1$.
• Figure 1: Function $\varepsilon \mapsto m(\varepsilon)$ for Example \ref{['ex:different_tau']}. We indicate as a blue dot the value of the function $m(\varepsilon^{\star})$, with $\varepsilon^\star = 0.7034$. We observe that Assumption \ref{['assumption_points']} holds, when $\varepsilon$ is sufficiently close to $\varepsilon^{\star}$.
• Figure 1: Contour plot for the smallest singular value of $\widehat{\mathcal{F}}_1(\lambda)+ \Delta \widehat{\mathcal{F}}_1(\lambda)$ at the grid of points obtained with meshgrid.
• Figure 2: Plot of the solution of the original problem (left) and perturbed problem (right) in the case of small delay.
• Figure 2: Additional check on the robustness of the approach, for Example \ref{['ex:different_tau']}, with $\tau=1$.

Theorems & Definitions (27)

• Remark 3.1
• Definition 3.2
• Remark 3.3
• Theorem 3.4
• Theorem 3.5
• Theorem 3.6
• Remark 4.1
• Definition 4.2
• Lemma 4.3: Constrained steepest descent method
• Proof 1