Structured distance to singularity as a nonlinear system of equations
March 05, 2026
Abstract
In this article we study the structured distance to singularity for a nonsingular matrix \(A\in\mathbb{C}^{n\times n}\), with a prescribed linear structure \(\mathcal{S}\) (for instance, a sparsity pattern, or a real Toeplitz structure), i.e., the norm of the smallest perturbation \(\Delta \in \mathcal{S}\), such that \(A + \Delta\) is singular. This is an example of structured matrix nearness problem: a family of problems that arise in control and systems theory and in numerical analysis, when characterizing the robustness of a certain property of a system with respect to perturbations that are constrained to a certain structure (for example the structure of the nominal system). We start by highlighting the parallelism between two main tools which have been proposed in the literature: a gradient system approach for a functional in the eigenvalues, which requires the solution of certain low-rank matrix differential equations (see [Guglielmi, Lubich, Sicilia, SINUM 2023]), and a two-level optimization approach in which the inner linear least-squares problem is solved explicitly (see [Usevich, Markovsky, JCAM 2014] and [Gnazzo, Noferini, Nyman, Poloni, FoCM 2025]). In particular, these articles underline the remarkable property that \(\Delta\) is (at least generically) the orthogonal projection onto the structure \(\mathcal{S}\) of a rank-1 matrix \(uv^*\). This property and the parallelism suggest a new reformulation of the problem into a system of nonlinear equations in the two vector unknowns \(u,v \in\mathbb{C}^n\). We study this new formulation, and propose an algorithm to solve these nonlinear equations directly with the multivariate Newton’s method. We discuss how to avoid the singularity of such system of nonlinear equations, and how to ensure monotonic convergence. The resulting algorithm is faster than the existing ones for large matrices, and maintains comparable accuracy.
Keywords: Structured matrix nearness problems, rank-1 perturbations, structured distance to singularity, nonlinear equations
AMS subject classifications: 15A18, 15A99, 65F15
1 Introduction↩︎
The aim of this work is to describe a new method to compute the nonsingularity radius of a structured matrix, which is a classical matrix nearness problem; for an extensive description of relevant nearness problems, we refer the reader to the seminal paper by N. Higham [1]. For a given nonsingular matrix \(A\), we wish to determine the size of the minimal additive perturbation which results in the loss of the considered property, that is, it makes it singular. In the unstructured case such a distance is equal to the minimal singular value of \(A\), as is stated by the well-known Eckart-Young-Mirsky theorem. However, when the matrix \(A\) has a certain structure, for example it is sparse or Toeplitz, which are cases where the structure is a linear manifold \(\mathcal{S}\), it appears more interesting to determine the closest singular matrix to \(A\) within the manifold \(\mathcal{S}\). In this case the Eckart-Young-Mirsky theorem is useless, since it provides generically an unstructured perturbation, and the mathematical (as well as the computational) problem is quite challenging. Our aim is to provide ideas and computational techniques to deal with this relevant problem, by proposing a novel algorithm merging ideas from [2] and [3].
For a general nonsingular complex matrix \(A \in \mathbb{C}^{n \times n}\), we consider the following problem: find a perturbation \(\Delta \in \mathbb{C}^{n \times n}\) of minimal norm such that \(A+\Delta\) becomes singular and \(\Delta \in \mathcal{S}\), where \(\mathcal{S}\) is a linear subspace in \(\mathbb{C}^{n\times n}\). The minimal norm of such a perturbation is called the nonsingularity radius of \(A\) \[\label{eq:distance} \min_{\Delta \in \mathcal{S}} \| \Delta \|^2 \ \text{s.t.}\ (A+\Delta) \text{ is singular}.\tag{1}\] In this article, the norm \(\|\cdot \|\) we consider in 1 is the Frobenius norm \(\|\cdot \|_F\) i.e., the Euclidean norm of the vector of the matrix entries: for \(M \in \mathbb{C}^{n \times n}\), its Frobenius norm is given by \[\| M \|_F = \Biggl( \sum_{i,j =1}^n |m_{ij}|^2 \Biggr)^{1/2}.\]
If we do not wish to preserve any structure, the nearest singular matrix to \(A\) is obtained directly from its singular value decomposition \[A = \sum_{i=1}^n \sigma_i\, u_i v_i^*,\] where \(u_i, v_i \in \mathbb{C}^n\) are the left and right singular vectors associated with the singular value \(\sigma_i > 0\), and \({}^*\) denotes the conjugate transpose. A nearest singular matrix is obtained by truncating the last term of the SVD expansion, i.e., by taking \[\Delta = - \sigma_n\, u_n v_n^*.\] Since the perturbation \(\Delta\) is a rank-one matrix, it is a minimizer also for the (operator) \(2\)-norm. Thus, the distance to singularity equals the smallest singular value \(\sigma_n\), for both the 2-norm and the Frobenius norm. This is a special case (rank \(r = n-1\)) of the classical Eckart–Young theorem [4], originally due to Schmidt [5]. However, when \(A\) belongs to a specified linear structure — for example, a fixed sparsity pattern, or the class of Toeplitz or Hankel matrices — the situation changes fundamentally: a nearest singular matrix within the same structure cannot be directly recovered by truncating the SVD of \(A\). A minimal structured perturbation \(\Delta\) is no longer rank one in general. Nevertheless, for the Frobenius norm we know that a minimal perturbation can be obtained as the orthogonal projection of a rank-one matrix onto the structure (see, e.g. [6], [3]).
This observation motivates the two-level iterative algorithm presented e.g. in [6], [3]. In the inner iteration we solve a system of differential equations for two vectors that depend on a distance parameter; this requires only matrix–vector multiplications with structured matrices and vector inner products. In the outer iteration we solve a scalar nonlinear equation to determine the structured distance to singularity. The gradient system approach has been proposed also in different contexts where matrix nearness problems arise: in stabilization theory [7], for computing the distance of coprime polynomials to common divisibility [8], in neural networks [9], in matrix polynomials [10] and in graph theory [11], [12], to name a few.
A different approach to solve the problem comes from the observation that if we impose the stricter constraint that \((A+\Delta)v=0\) for a prescribed vector \(v\), the optimal \(\Delta_*(v)\) can be found explicitly. Hence, one need only to optimize on \(v\). This approach has been studied extensively in a series of papers by Usevich and Markovsky [13], [14], even though this idea of “variable projection” dates back at least to [15]. The approach has recently been revisited and extended in [2], in particular incorporating a regularization technique. Also this technique can be applied to many of the problems described above, such as stabilization, matrix polynomials, and polynomial GCD.
Outline↩︎
The paper is organized as follows. In Section 2 we briefly describe two existing methods for the numerical approximation of the structured distance to singularity, namely the ODE–based method and the structured low-rank approximation (SLRA)/Riemann-Oracle method. Moreover, we provide a list of parallels between the two optimization procedures. In Section 3, we introduce a novel approach for the computation of the structured distance to singularity, recasting the problem as a system of nonlinear equations. In Section 4, we describe the implementation details, and in Section 5, we test the behavior of the proposed method via numerical experiments. Finally, Section 6 summarizes the main results.
2 The two existing methods↩︎
2.1 A gradient system approach↩︎
In this paragraph we describe the gradient system approach proposed in many recent works (see e.g. [3], [16]–[18]) for solving matrix nearness problems. We consider its general framework and we focus on how it is applied to the problem of computing the structured distance to singularity.
Given a matrix \(A\in \mathbb{C}^{n\times n}\) and a property \(\mathscr{P}\) related to the spectrum of \(A\) we look for the solution of the optimization problem \[\label{prob:dist95stru} \mathop{\mathrm{arg\,min}}_{\Delta \in \mathcal{S}}\{\|\Delta\|_F : A+\Delta \textrm{ does not fulfil the property } \mathscr{P}\},\tag{2}\] where \(\mathcal{S}\subseteq \mathbb{C}^{n\times n}\) is a subspace that enforces a structure on the sought perturbation \(\Delta\). The gradient system approach relies on a two-level method that splits the original problem into two nested sub-problems that are solved by an inner iteration and by an outer iteration. In the general setting, we consider a functional \(\mathscr{F}:\mathcal{S}\rightarrow \mathbb{R}\) such that \(\mathscr{F}(0)>0\) and, for all \(\Delta\in\mathcal{S}\), \[\mathscr{F}(\Delta)\leq 0 \iff A+\Delta \textrm{ does not fulfil the property } \mathscr{P},\] which means that the functional \(\mathscr{F}\) takes non-positive values if and only if \(\Delta\) is admissible. We rewrite the perturbation \(\Delta=\varepsilon E\), where \(\varepsilon>0\) is the perturbation size and \(E\) has unit Frobenius norm and we define the functional \(F_\varepsilon\) as \[F_\varepsilon(E):=\mathscr{F}(\varepsilon E).\] The inner iteration minimizes the functional \(F_\varepsilon\) when the perturbation size \(\varepsilon\) is fixed, while the outer iteration aims to find the smallest value \(\varepsilon_\star\) such that it is possible to have \(\mathscr{F}(\Delta)= 0\) for some perturbation \(\Delta\in \mathcal{S}\) with Frobenius norm \(\varepsilon_\star\). The outline of the two-level method is the following:
Inner iteration: For a fixed \(\varepsilon\), compute a matrix perturbation \(E_\star(\varepsilon)\) such that \[\label{prob:inn95iter95stru} E_\star(\varepsilon)\in\mathop{\mathrm{arg\,min}}_{\|E\|_F=1, \;E\in \mathcal{S}} F_\varepsilon(E)=\mathop{\mathrm{arg\,min}}_{\|\Delta\|_F=\varepsilon, \;\Delta \in \mathcal{S}} \mathscr{F}(\Delta),\tag{3}\]
Outer Iteration: Find the smallest value \(\varepsilon_\star>0\) such that \[\label{eq:zero} \phi(\varepsilon):=F_\varepsilon(E_\star(\varepsilon))=0.\tag{4}\]
The inner iteration is the most elaborated procedure, while the outer iteration is theoretically easier to solve, since it consists of a one-dimensional root-finding problem, even though it could still be challenging.
In this work we consider the approach of [3] used for computing the structured distance to singularity, meaning that the property we aim to violate is \(\mathscr{P}=\{\textrm{the matrix is nonsingular}\}\) and the functional \(\mathscr{F}\) takes the form \[\label{eigenvalue32optimization32problem} \mathscr{F}(\Delta)=|\lambda_{\textrm{target}}(A+\Delta)|^2,\tag{5}\] where \(\lambda_{\textrm{target}}\) is the eigenvalue with smallest absolute value. We also restrict to consider the case where \(\mathcal{S}\) describes the sparsity pattern of a matrix. To solve problem 3 , the inner iteration introduces a perturbation matrix path \(E(t)\) with \(t\geq 0\) and it integrates the matrix ODE \[\label{eq:odeEstru} \dot{E}=-\Pi_\mathcal{S}G_\varepsilon(E)+\mathop{\mathrm{\textrm{Re}}}\langle \Pi_\mathcal{S}G_\varepsilon(E),E\rangle E,\tag{6}\] where - for two matrices \(M,N\) - we denote by \[\scalar{M,N} = \operatorname{Tr}(M^*N)\] the scalar product on matrices that induces the Frobenius norm, by \(G_\varepsilon(E)\) the gradient of \(F_\varepsilon\), and by \(\Pi_\mathcal{S}\) the orthogonal projection onto \(\mathcal{S}\). The expression of \(G_\varepsilon(E)\) is \[\label{GepsExy} G_\varepsilon(E)=-\lambda xy^*,\tag{7}\] where \(x\) and \(y\) are, respectively, the unit left and right eigenvectors associated with the target eigenvalue \(\lambda=\lambda_{\textrm{target}}(A+\varepsilon E)\) so that \(x^*y>0\); the orthogonal projection \(\Pi_\mathcal{S}\) simply consists in replacing by \(0\) the entries outside of the pattern.
It is possible to prove that the stationary points of 6 corresponds to the local minima of \(F_\varepsilon\) and, in order to find them, we integrate the ODE 6 until we reach a sought stationary point. Since equation 6 is a gradient system, the integration will always end up in a stationary point and, up to non-generic events, these have the form \(E\propto \Pi_\mathcal{S}G_\varepsilon(E)\).
The matrix \(G_\varepsilon(E)\) has rank-\(1\) and hence it follows that the stationary points are the projections onto \(\mathcal{S}\) of a rank-\(1\) matrix. This motivates to consider a different ODE to solve the inner iteration, whose trajectory belongs to the rank-\(1\) manifold \(\mathcal{M}_1\): introduce a rank-\(1\) matrix path \(Y(t)\subseteq \mathcal{M}_1\) such that \[E(t)=\Pi_\mathcal{S}Y(t), \qquad t\in [0,+\infty)\] and we consider the ODE \[\label{eq:odeY} \dot{Y}=P_Y\left(-G_\varepsilon(\Pi_\mathcal{S}Y)+\mathop{\mathrm{\textrm{Re}}}\langle P_Y(G_\varepsilon(\Pi_\mathcal{S}Y)),\Pi_\mathcal{S}Y \rangle Y\right),\tag{8}\] where \(P_Y\) denotes the orthogonal projection with respect to the Frobenius inner product onto the tangent space \(\mathcal{T}_Y \mathcal{M}_1\) of \(\mathcal{M}_1\) in \(Y\). There exists an explicit one-to-one correspondence between the stationary points of 6 and those of 8 (see [3]), which ensures that we are not introducing nor losing solutions of problem 2 if we integrate the low-rank equation instead of the full-rank one. The rank-1 differential equation 8 for \(Y=\rho \hat{u} v^*\) can be restated in terms of differential equations for the unit norm vectors \(\hat{u}, v\) and an explicit formula for \(\rho\).
Lemma 1 (Differential equations for the factors). Every solution \(Y(t)\in \mathcal{M}_1\) of the rank-1 differential equation 8 with \(\| \Pi_\mathcal{S}Y(t) \|_F=1\) can be written as \(Y(t)=\rho(t)\hat{u}(t)v(t)^*\) where \(\hat{u}(t)\) and \(v(t)\) of unit norm satisfy the differential equations \[\begin{align} \rho \dot{\hat{u}} &= - \displaystyle \tfrac\mathrm{i}2 \mathop{\mathrm{\textrm{Im}}}(\hat{u}^*Gv)\hat{u} -(I-\hat{u}\hat{u}^*) Gv , \qquad \\ \rho \dot{v} &= - \displaystyle \tfrac\mathrm{i}2 \mathop{\mathrm{\textrm{Im}}}(v^*G\hat{u})v -(I-vv^*) G^*\hat{u} , \end{align}\] where \(G=G_\varepsilon(E)\) for \(E=\Pi_\mathcal{S}Y = \rho\, \Pi_\mathcal{S}(\hat{u}v^*)\) and \(\rho=1/\| \Pi_\mathcal{S}(\hat{u}v^*)\|_F\).
Even though equation 8 is not a gradient system, it is somehow close to being one: see [3]. In particular, when choosing a proper starting point sufficiently close to a stationary point, integrating equation 8 always leads to that stationary point. This just yields a local convergence result, weaker than the global convergence property of a gradient system.
The first observation is that at a stationary point \(Y\) of 8 , we have \(P_Y G_\varepsilon(E)=G_\varepsilon(E)\) for \(E=\Pi_\mathcal{S}Y\). Therefore, close to a stationary point, \(P_Y G_\varepsilon(E)\) will be close to \(G_\varepsilon(E)\). It turns out that it is even quadratically close, as is stated in the following lemma.
Lemma 2 (Projected gradient near a stationary point).
Let \(Y_\star\in \mathcal{M}_1\) with \(E_\star=\Pi_\mathcal{S}Y_\star \in \mathcal{S}\) of unit Frobenius norm. Let \(Y_\star\) be a stationary point of the rank-1 projected differential equation 8 , with an associated target eigenvalue \(\lambda\) of \(A+\varepsilon E_\star\) that is simple. Then, there exist \(\bar \delta>0\) and a real \(C\) such that for all positive \(\delta\le \bar\delta\) and all \(Y \in \mathcal{M}_1\) with \(\| Y-Y_\star \| \le \delta\) and associated \(E=\Pi_\mathcal{S}Y\) of unit norm, we have \[\| P_Y G_\varepsilon\left(E \right) - G_\varepsilon\left(E \right) \| \le C \delta^2.\]
As a direct consequence of this lemma, a comparison of the differential equations 8 and 6 yields that when \(\delta\)-close to a stationary point, the functional decreases monotonically along solutions of 8 up to \(O(\delta^2)\), and even with the same negative derivative as for the gradient flow 6 up to \(O(\delta^2)\). Note that the derivative of the functional is proportional to \(-\delta\) in a \(\delta\)-neighbourhood of a strong local minimum. Guglielmi, Lubich & Sicilia [3] used Lemma 2 to prove a result on local convergence as \(t\to\infty\) to strong local minima of the functional \(F_\varepsilon\) for \(E(t)=\Pi_\mathcal{S}Y(t)\) of unit Frobenius norm associated with solutions \(Y(t)\) of the rank-1 differential equation 8 .
Numerical experiments show that this is enough to make the method work in practice. In this way, integrating 8 makes it possible to exploit the underlying low-rank features of the matrix nearness problem getting some benefits in the numerical computations.
Remark 1. In some frameworks, the eigenvalue optimization problem 5 is replaced by the singular value optimization problem \[\widetilde{\mathscr{F}}(\Delta)=\sigma_{\textrm{target}}(A+\Delta),\] where we substitute \(\lambda_\textrm{target}^2\) by \(\sigma_\textrm{target}\) in the definition of \(\mathscr{F}\). This idea is similar to the approach previously employed in [19], for instance. In this paper we will consider this specific case, and in particular we will focus on the smallest singular value, that is \(\sigma_\textrm{target}=\sigma_{\min}\), so that \[\mathop{\mathrm{arg\,min}}_{\Delta\in \mathcal{S}} \widetilde{\mathscr{F}}(A+\Delta)\] provides the structured distance to singularity of \(A\). In this case it is also possible to study the corresponding gradient system similar to 6 , i.e. \[\label{eq:odeEstru2} \dot{E}=-\Pi_\mathcal{S}\widetilde{G}_\varepsilon(E)+\mathop{\mathrm{\textrm{Re}}}\langle \Pi_\mathcal{S}\widetilde{G}_\varepsilon(E),E\rangle\, E,\tag{9}\] where \(\widetilde{G}=u v^*\) is the gradient associated with the function \(\widetilde{F}_\varepsilon(E):=\widetilde{\mathscr{F}}(\varepsilon E)\) (analogous to \(F_\varepsilon(E)\) for the eigenvalue case) and \(u\) and \(v\) are the left and right singular vectors associated with \(\sigma_{\min}\)
2.2 SLRA / Riemann-Oracle approach↩︎
In this section, we describe a different framework that was studied in a series of papers by Markovsky and Usevich [13], [14], [20], with the name of structured low-rank approximation, and then expanded to more general nearness problems and studied in more detail in [2], with the name of Riemann oracle. For the purpose of this work, we briefly describe the method for the specific case of the computation of the structured distance to singularity for a given nonsingular matrix. We consider again a matrix \(A\in \mathbb{C}^{n\times n}\) and a linear subspace \(\mathcal{S}\) of dimension \(p\) containing the admissible perturbations \(\Delta\in \mathbb{C}^{n\times n}\). Let \(P^{(1)}, \dots, P^{(p)} \in \mathbb{C}^{n\times n}\) be an orthonormal basis of \(\mathcal{S}\), i.e., \[\mathcal{S} = \left\{ \Delta \in \mathbb{C}^{n\times n} \colon \Delta = \sum_{i=1}^p P^{(i)} \delta_i, \quad \delta = \begin{bmatrix} \delta_1\\ \vdots\\ \delta_p \end{bmatrix} \in \mathbb{C}^p \right\}, \quad \operatorname{Tr}((P^{(i)})^*P^{(j)}) = \delta_{ij}.\] Equivalently, \(\mathop{\mathrm{vec}}(\Delta) = \mathcal{P}\delta\), where \[\label{eq:matrix95P} \mathcal{P} = \begin{bmatrix} \mathop{\mathrm{vec}}P^{(1)} & \dots & \mathop{\mathrm{vec}}P^{(p)} \end{bmatrix} \in \mathbb{C}^{n^2\times p}\tag{10}\] has orthonormal columns.
The method is based on the observation that the problem 1 becomes easier if we fix a vector \(v\in\mathbb{C}, v\neq 0\) that must be in the kernel of \(A+\Delta\) (in the formulation of [2], this target vector is provided by an oracle). To solve \[\label{minfv} \min \, \|\Delta\|_F^2 \ \text{s.t.}\ {\Delta \in \mathcal{S}}, \, (A+\Delta)v = 0,\tag{11}\] we can eliminate the constraint \(\Delta\in\mathcal{S}\) by rewriting 11 in terms of the vector \(\delta\) such that \(\mathop{\mathrm{vec}}\Delta = \mathcal{P}\delta\); then 11 becomes \[\label{minfvdelta} \min\, \norm{\delta}^2 \ \text{s.t.}\ M\delta = r,\tag{12}\] where \(r = r(v) = -Av\in\mathbb{C}^n\) and \[\label{Mr} M = M(v) = (v^\top\otimes I_n)\mathcal{P} = \begin{bmatrix} P^{(1)}v & P^{(2)}v & \dots & P^{(p)}v \end{bmatrix} \in \mathbb{C}^{n \times p}.\tag{13}\] This is a classical problem whose (unique) solution is \(\delta_* = M^+r\), where \(M^+\) denotes the Moore-Penrose pseudoinverse (assuming that the feasible region is non-empty, i.e., \(r\in \text{range}(M)\)).
Hence, we know how to solve the inner minimization subproblem over \(\Delta\) in closed form and compute the optimal \(\norm{\Delta}_F\) for a given candidate kernel vector \(v\), which we can take in the unit sphere in \(\mathbb{C}^n\), here denoted by \(\mathbb{S}_1\). So to solve the original problem 1 for the Frobenius norm, we solve the outer minimization problem \[\min_{v \in \mathbb{S}_1} f(v), \quad f(v) = \norm{M(v)^+r(v)}^2,\] via Riemannian optimization over the manifold \(\mathbb{S}_1\).
In [14], these formulas are obtained under the additional assumption that \(M\) has full column rank. However, as noted in [2], this assumption is problematic: when the rank of \(M(v)\) drops, the function \(f(v)\) has a removable discontinuity; and the global minimum of \(f(v)\) may occur in this discontinuity point. In this case, a classical derivative-based optimization method for \(f(v)\) would miss this discontinuity and return instead a suboptimal local minimum.
In [2], this issue is solved using a regularization procedure: we minimize a relaxed version of the objective functional \(f(v)\), namely \[\label{fepsilonv} f_\varepsilon(v) = \min_{{\Delta\in\mathcal{S}}} \,\|\Delta\|_F^2 + \varepsilon^{-1} \|(A+\Delta)v\|^2,\tag{14}\] for a given \(\varepsilon>0\). Arguing as above, one gets closed-form expressions for the minimum \[f_\varepsilon(v) = r^*(MM^* + \varepsilon I)^{-1}r\] and for the corresponding argument minimum \[\label{deltastar} \delta_* = M^*(MM^*+\varepsilon I)^{-1}r = (M^*M+\varepsilon I)^{-1}M^*r, \quad \mathop{\mathrm{vec}}\Delta_* = \mathcal{P}\delta_*.\tag{15}\] When \(\varepsilon\to 0\), the minimum of \(f_\varepsilon(v)\) tends to the minimum of \(f(v)\) [2].
Moreover, following the approach used in [21] for a similar problem, we observe that \(\mathcal{P}\mathcal{P}^*\) is the orthogonal projection matrix over \(\mathop{\mathrm{vec}}\mathcal{S}\), hence if we set \(u = (MM^*+\varepsilon I)^{-1}r\in\mathbb{C}^n\), we have \(\delta_* = M^*u\) and \[\label{vecDelta} \mathop{\mathrm{vec}}(\Delta_*) = \mathcal{P}\delta_* = \mathcal{P}M^*u = \mathcal{P}\mathcal{P}^*(\overline{v} \otimes u) = \mathop{\mathrm{vec}}\Pi_\mathcal{S}(uv^*),\tag{16}\] i.e., the matrix \(\Delta_*\) computed at each iteration is the projection on the structure \(\mathcal{S}\) of the rank-1 matrix \(uv^*\).
The Euclidean gradient of 14 is \[\nabla_v f_\varepsilon= (A+\Delta)^*u, \quad \Delta = \Pi_\mathcal{S}(uv^*).\] In a stationary point \(v_*\) of \(f_\varepsilon\) on the unit sphere, the Euclidean gradient \(\nabla_{v_*} f\) must be a multiple of \(v_*\). In the limit \(\varepsilon\to 0\) the function \(f(v)\) is scale-invariant, i.e., \(f(v\alpha)=f(v)\) for each \(\alpha\in\mathbb{C},v\in\mathbb{C}^n\), so the radial component of the gradient vanishes, and it must be the case that \(\nabla_{v_*} f = 0\).
As mentioned, in this work we are mainly interested in sparsity structures. In this case, \(MM^*\) becomes a diagonal matrix (see [2] and the proof of Lemma 5 below), making the method simpler and faster.
2.3 Similarities between the two methods↩︎
It is interesting to note the parallels between the two methods.
| ODE approach | Riemann–Oracle approach |
|---|---|
| Inner iteration: for given \(\varepsilon>0\), use an ODE integrator to compute \(\Delta_* = \Pi_\mathcal{S}(uv^*)\), with \(u=\varepsilon\rho \hat{u}\), which satisfies \(\dot{Y}=0\). | Inner iteration: for given \(\varepsilon> 0\), use Riemannian optimization to compute \(\Delta_* = \Pi_{\mathcal{S}}(uv^*)\) which satisfies \(\nabla_v f_{\varepsilon} = v\mu\). |
| Outer iteration: compute optimal \(\varepsilon\) by Newton-bisection to obtain \(\lambda = 0\) (univariate minimization problem). | Outer iteration: let \(\varepsilon\to 0\). |
| At each inner iteration: \((A+\Delta_*)v=v\lambda\), \((A+\Delta_*)^*u = u \overline{\lambda}\). | At each inner iteration: \(\|\Delta\|_F^2 + \varepsilon^{-1} \|(A+\Delta)v\|^2\) is minimized and \((A+\Delta_*)^*u = v\mu\). |
| At convergence of the outer iteration: \((A+\Delta_*)v=0\), \(u^*(A+\Delta_*) = 0\). | |
While the two methods are not the same, their structure is very similar: they minimize different functions at each step, but in the end of their nested iterations they compute a pair \((u,v)\) which satisfies the same nonlinear system of equations.
3 Reformulation as a nonlinear system↩︎
These two formulations suggest looking for vectors \(u,v\in\mathbb{C}^n\) that satisfy the two required equations directly.
Problem: find \(u,v\in\mathbb{C}^n\) that satisfy the non-linear system of equations \[\label{nlsys} (A + \Pi_\mathcal{S}(uv^*))v = (A + \Pi_\mathcal{S}(uv^*))^* u = 0.\tag{17}\]
When these two equations 17 hold, the pair \((u,v)\) is a critical point of the optimization problems in both the ODE and the Oracle approach.
We prove the following statement: the vectors \(u,v\) computed by the ODE approach are, up to a rescaling, a solution of Problem 17 .
Theorem 1 (Stationary points for fixed \(\varepsilon\)). Let \(\widetilde{\mathscr{F}}(\Delta)=\sigma_{\min}(A+\Delta)\) and fix \(\varepsilon>0\). Assume that \(\sigma_{\min}(A+\varepsilon E)\) is positive and simple along the trajectory considered by the ODE approach. Then any stationary point \(E(\varepsilon)\in \mathcal{S}\) of the constrained flow with \(\|E(\varepsilon)\|_F=1\) satisfies, up to multiplication by a real scalar, \[E(\varepsilon) \propto \Pi_\mathcal{S}\!\bigl(u(\varepsilon)v(\varepsilon)^*\bigr), \label{eq:rank1proj-stationary}\tag{18}\] where \(u(\varepsilon),v(\varepsilon)\) are the left and right singular vectors of unit norm, associated with \(\sigma(\varepsilon):=\sigma_{\min}(A+\varepsilon E(\varepsilon))\). In particular, with the normalization \[E(\varepsilon) = \frac{\Pi_\mathcal{S}(u(\varepsilon)v(\varepsilon)^*)}{\|\Pi_\mathcal{S}(u(\varepsilon)v(\varepsilon)^*)\|_F}, \label{eq:E-structure-thm1}\tag{19}\] the singular-vector relations read \[(A+\varepsilon E(\varepsilon))v(\varepsilon)=\sigma(\varepsilon)u(\varepsilon), \qquad (A+\varepsilon E(\varepsilon))^*u(\varepsilon)=\sigma(\varepsilon)v(\varepsilon). \label{eq:sv-rel-thm1}\tag{20}\]
Proof. Let \(\Delta(t)=\varepsilon E(t)\) with \(\|E(t)\|_F\equiv 1\). As in Remark 1, we define the functional \(\widetilde{F}_\varepsilon(E)=\widetilde{\mathscr{F}}(\varepsilon E)=\sigma_{\min}(A+\varepsilon E)\). Whenever \(\sigma_{\min}\) is simple and strictly positive one has the standard differential identity \[\frac{\mathrm{d}}{\mathrm{d}t}\widetilde{F}_\varepsilon(E(t)) = \varepsilon\,\mathop{\mathrm{\textrm{Re}}}\langle \dot{E}(t),\,u(t)v(t)^*\rangle,\] where \(u(t),v(t)\) are the corresponding left/right singular vectors. Hence the (Euclidean) gradient of \(\widetilde{F}_\varepsilon\) at \(E\) is \(\widetilde{G}_\varepsilon(E)=u v^*\).
The constrained gradient flow on the Frobenius-norm unit sphere in \(S\) is \[\dot{E}=-\Pi_\mathcal{S}\widetilde{G}_\varepsilon(E)+\mathop{\mathrm{\textrm{Re}}}\langle \Pi_\mathcal{S}\widetilde{G}_\varepsilon(E),E\rangle\,E,\] so at a stationary point \(E(\varepsilon)\) we have \(-\Pi_\mathcal{S}\widetilde{G}_\varepsilon(E(\varepsilon))+\mathop{\mathrm{\textrm{Re}}}\langle \Pi_\mathcal{S}\widetilde{G}_\varepsilon(E(\varepsilon)), E(\varepsilon)\rangle E(\varepsilon)=0\), which implies 18 . Choosing the normalization 19 yields 20 . ◻
Note that we had to assume \(\sigma_{\min}>0\) in the previous theorem. In the next theorem we extend to the limit to obtain a result for \(\sigma_{\min}=0\).
Theorem 2 (Limit characterization at \(\varepsilon_\star\)). Assume that for every \(\varepsilon \in (0,\varepsilon_\star)\) there exists a stationary point \(E(\varepsilon)\in \mathcal{S}\) of unit Frobenius norm, \(\|E(\varepsilon)\|_F = 1\), such that the smallest singular value \[\sigma(\varepsilon) := \sigma_{\min}(A+\varepsilon E(\varepsilon))\] is positive and simple, and the associated left and right singular vectors \(u(\varepsilon), v(\varepsilon)\) (of unit norm) satisfy \[(A+\varepsilon E(\varepsilon)) v(\varepsilon) = \sigma(\varepsilon) u(\varepsilon), \qquad (A+\varepsilon E(\varepsilon))^* u(\varepsilon) = \sigma(\varepsilon) v(\varepsilon), \label{eq:sv-rel}\tag{21}\] and \[E(\varepsilon) = \frac{\Pi_\mathcal{S}(u(\varepsilon)v(\varepsilon)^*)}{\|\Pi_\mathcal{S}(u(\varepsilon)v(\varepsilon)^*)\|_F}. \label{eq:E-structure}\tag{22}\] Moreover suppose that \[\lim_{\varepsilon \nearrow \varepsilon_\star} \sigma(\varepsilon) = 0,\] and that there exists a constant \(c>0\) such that \[\|\Pi_\mathcal{S}(u(\varepsilon)v(\varepsilon)^*)\|_F \ge c \qquad \text{for all } \varepsilon\in(0,\varepsilon_\star). \label{eq:nondeg}\tag{23}\] Then there exist vectors \(\widehat{u}, \widehat{v} \in \mathbb{C}^n\), with \(\|\widehat{v}\|_2=1\), and a matrix \(\widehat{\Delta} := \Pi_\mathcal{S}(\widehat{u}{\widehat{v}}^*)\) such that \[(A+\widehat{\Delta}) \widehat{v} = 0, \qquad and \qquad (A+\widehat{\Delta})^* \widehat{u} = 0, \label{eq:limit-system}\tag{24}\] i.e., \((\widehat{u},\widehat{v})\) satisfies the nonlinear system \((A+\Pi_\mathcal{S}(\widehat{u}{\widehat{v}}^*))\widehat{v} = (A+\Pi_\mathcal{S}(\widehat{u}{\widehat{v}}^*))^*\widehat{u} = 0.\)
Proof. Define \(\Delta(\varepsilon) := \varepsilon E(\varepsilon)\). By 22 we can rewrite \[\Delta(\varepsilon) = \Pi_\mathcal{S}\!\left( \alpha(\varepsilon)\, u(\varepsilon)v(\varepsilon)^* \right), \qquad where \qquad \alpha(\varepsilon) := \frac{\varepsilon}{\|\Pi_\mathcal{S}(u(\varepsilon)v(\varepsilon)^*)\|_F}.\] Set \(\tilde{u}(\varepsilon) := \alpha(\varepsilon)\, u(\varepsilon)\). Then \[\Delta(\varepsilon) = \Pi_\mathcal{S}(\tilde{u}(\varepsilon)v(\varepsilon)^*). \label{eq:Delta-rewrite}\tag{25}\] By 23 , the scalars \(\alpha(\varepsilon)\) are bounded on \((0,\varepsilon_\star)\), hence \(\tilde{u}(\varepsilon)\) is bounded as well, since \(\|u(\varepsilon)\|_2=1\). Moreover, \(\|v(\varepsilon)\|_2=1\) for all \(\varepsilon\). Therefore, by compactness of the unit sphere, there exists a sequence \(\varepsilon_k \nearrow \varepsilon_\star\) such that \[v(\varepsilon_k) \to \widehat{v}, \qquad \tilde{u}(\varepsilon_k) \to \widehat{u}\] for some \(\widehat{u}, \widehat{v}\) with \(\|\widehat{v}\|_2=1\). Since \(\|\Delta(\varepsilon)\|_F=\varepsilon\), the sequence \(\Delta(\varepsilon_k)\) is bounded and hence (up to subsequences) \(\Delta(\varepsilon_k) \to \widehat{\Delta} \in \mathcal{S}.\) Passing to the limit in 25 and using continuity of \(\Pi_\mathcal{S}\) yields \[\label{eq:Deltahat} \widehat{\Delta} = \Pi_\mathcal{S}(\widehat{u}{\widehat{v}}^*).\tag{26}\]
Now consider the first relation in 21 : \((A+\Delta(\varepsilon_k)) v(\varepsilon_k) = \sigma(\varepsilon_k) u(\varepsilon_k)\). Taking norms gives \[\|(A+\Delta(\varepsilon_k)) v(\varepsilon_k)\|_2 = \sigma(\varepsilon_k),\] and by assumption \(\sigma(\varepsilon_k)\to 0\). Since \[A+\Delta(\varepsilon_k) \to A+\widehat{\Delta}, \qquad v(\varepsilon_k) \to \widehat{v},\] we obtain \((A+\widehat{\Delta}) \widehat{v} = 0\).
For the adjoint relation, multiply the second equation in 21 by \(\alpha(\varepsilon_k)\): \[(A+\Delta(\varepsilon_k))^* \tilde{u}(\varepsilon_k) = \alpha(\varepsilon_k)\sigma(\varepsilon_k) \, v(\varepsilon_k).\] The right-hand side tends to zero because \(\alpha(\varepsilon_k)\) is bounded and \(\sigma(\varepsilon_k)\to 0\). Passing to the limit yields \((A+\widehat{\Delta})^* \widehat{u} = 0\). This proves 24 . ◻
The proof can be simplified if we assume that \(\exists \lim\limits_{\varepsilon\nearrow \varepsilon_\star} u(\varepsilon) = \widehat{u}\) and \(\lim\limits_{\varepsilon\nearrow \varepsilon_\star} v(\varepsilon) = \widehat{v}\).
We can also prove an explicit stationarity property of the limit matrix.
Theorem 3 (Clarke-stationarity at \(\widehat\Delta\)). Under the assumptions of Theorem 2, let \[\mathcal{S}_{\varepsilon_\star}:=\{\Delta\in \mathcal{S}:\;\|\Delta\|_F=\varepsilon_\star\}, \qquad T_{\widehat\Delta}\mathcal{S}_{\varepsilon_\star} =\{Z\in \mathcal{S}:\;\mathop{\mathrm{\textrm{Re}}}\langle Z,\widehat\Delta\rangle=0\}\] and \[\widetilde{\mathscr{F}}(\Delta):=\sigma_{\min}(A+\Delta),\qquad \Delta\in \mathcal{S}.\]
Then \(\widehat \Delta\) is Clarke stationary (see [22] for more details) on \(\mathcal{S}_{\varepsilon_\star}\). Equivalently, with \(\Pi\) the projector onto the tangent space, \[0\in \Pi_{T_{\widehat\Delta}\mathcal{S}_{\varepsilon_\star}} \bigl(\partial^C \widetilde{\mathscr{F}}(\widehat\Delta)\bigr),\] i.e.there exists \(G\in \partial^C\sigma_{\min}(A+\widehat\Delta)\) (the Clarke subgradient) such that \(\Pi_\mathcal{S}(G)\) is collinear with \(\widehat\Delta\).
Proof. Using the same arguments of Theorem 2, we let \(\Delta_k = \varepsilon_k E_k\). Since \(\sigma_k = \sigma_{\min}(A + \Delta_k) >0\) is simple, \(X\mapsto \sigma_{\min}(X)\) is differentiable at \(X_k:=A+\Delta_k\) and its gradient is \(u_k v_k^*\). Stationarity of \(E_k\) on the unit sphere implies the gradient is collinear with \(E_k\).
Since \(\|u_k v_k^*\|_F=1\), the sequence \(\{u_k v_k^*\}\) is bounded; consider a convergent subsequence such that \(u_k v_k^*\to G\). We obtain \[\Pi_\mathcal{S}(u_k v_k^*)=\|\Pi_\mathcal{S}(u_k v_k^*)\|_F\,E_k \;\Longrightarrow\; \Pi_\mathcal{S}(G)=\eta\,\widehat E, \qquad \widehat E:=\widehat\Delta/\varepsilon_\star,\] for some \(\eta > 0\) (by the nondegeneracy assumption 23 ), hence \(\Pi_\mathcal{S}(G)\) is collinear with \(\widehat\Delta\). The map \(X\mapsto \sigma_{\min}(X)\) is locally Lipschitz, and its Clarke subdifferential is outer semicontinuous. Thus, since \(X_k\to \widehat X:=A+\widehat\Delta\) and \(u_k v_k^*=\nabla\sigma_{\min}(X_k)\), any limit \(G\) of \(u_k v_k^*\) belongs to \(\partial^C\sigma_{\min}(\widehat X)\). Therefore \(G\in\partial^C \widetilde{\mathscr{F}}(\widehat\Delta)\) and \(\Pi_\mathcal{S}(G)\) is collinear with \(\widehat\Delta\), proving the Clarke stationarity statement. ◻
Remark 2. If \((u,v)\) is a solution to 17 , then \[\label{normalization32ambiguity} (u\alpha^{-1}, v\overline{\alpha})\tag{27}\] is another solution for each \(\alpha\in\mathbb{C}, \alpha\neq 0\); hence the solutions are defined up to a normalization factor. We can assume without loss of generality that \(\norm{v}=1\).
3.1 Differentials, gradients and Hessians↩︎
We set \[G(u,v) = \begin{bmatrix} (A+\Delta)v\\ (A+\Delta)^*u \end{bmatrix}, \quad \Delta = \Pi_\mathcal{S}(uv^*),\] so that 17 becomes \(G(u,v) = 0\). It is interesting to note that, when \(A\in\mathcal{S}\), this quantity \(G(u,v)\) is the gradient of a scalar function \(F(u,v)\). We prove it in the following lemma.
Lemma 3. Consider the function \(F(u,v) = \frac{1}{2}\norm{A+\Pi_\mathcal{S}(uv^*)}_F^2\). Then, the two partial gradients of \(F\) with respect to \(u\) and \(v\) are \[\label{generalgradient} \nabla_u F(u,v) = (\Pi_\mathcal{S}(A)+\Delta)v, \quad \nabla_v F(u,v) = (\Pi_\mathcal{S}(A)+\Delta)^*u, \quad \Delta = \Pi_\mathcal{S}(uv^*).\tag{28}\]
Proof. Note that \[\scalar{A+\Delta, A+\Delta} = \|A\|_F^2 + 2 \scalar{\Pi_\mathcal{S}(A),uv^*}+\|\Pi_\mathcal{S}(uv^*)\|_F^2.\] We compute the differential \[\begin{align} dF(u,v) &= d \frac{1}{2}\left(\|A\|_F^2 + 2 \scalar{\Pi_\mathcal{S}(A),uv^*}+\|\Pi_\mathcal{S}(uv^*)\|_F^2\right)\\ &= \scalar{\Pi_\mathcal{S}(A), (du) v^* + u (dv)^*}+\scalar{\Pi_\mathcal{S}((du) v^* + u (dv)^*),\Pi_\mathcal{S}(uv^*)}\\ &= \scalar{\Pi_\mathcal{S}(A)+\Delta, (du) v^* + u (dv)^*}\\ &= \scalar{(\Pi_\mathcal{S}(A)+\Delta)v, du} + \scalar{dv, (\Pi_\mathcal{S}(A)+\Delta)^*u}, \end{align}\] where we have used the definition of \(\scalar{\cdot,\cdot}\), the fact that \(\Pi_\mathcal{S}\) is an orthogonal projection, and the cyclic property of the trace. From the last line we can read off the two gradients. ◻
When \(A \in \mathcal{S}\), \(\Pi_\mathcal{S}(A) = A\) and hence 28 coincide with the two blocks of \(G(u,v)\). This result shows that 17 is a gradient system.
The relation between this \(F(u,v)\) and the original minimization problem 1 is not immediate; and in general, the global distance minimizer \(\Delta_*\) is neither a minimum nor a maximum of the functional \(F\): see Example 1 below.
We can give an explicit formula for the differential \(dG(u,v)\) of \(G(u,v)\), which is a \(\mathbb{R}\)-linear operator \(H: \mathbb{C}^{2n} \to \mathbb{C}^{2n}\). We use the letter \(H\) because, in the case in which \(G\) is the gradient of \(F(u,v)\), the operator \(H\) is the Hessian of \(F(u,v)\).
We note that, even when \(u,v\) are complex vectors, \(H\) is only guaranteed to be \(\mathbb{R}\)-linear, since conjugates will appear in its expression. To fix the notation, let \(\overline{x}\) denote the (entrywise) conjugate of a scalar, vector or matrix \(x\); while \(x^\top\) stands for the transpose of \(x\).
Lemma 4. The differential of \(G(u,v)\) is the \(\mathbb{R}\)-linear operator \(H\) such that for each \(\begin{bsmallmatrix} du\\dv \end{bsmallmatrix} \in \mathbb{C}^{2n}\) \[H \begin{bmatrix} du\\ dv \end{bmatrix} = \begin{bmatrix} MM^* \, du + (A+\Delta)dv + M N^\top\overline{dv}\\ (A+\Delta)^*du + NM^\top \overline{du} + NN^*\,dv \end{bmatrix},\] where \(\mathcal{P}\) is defined as in 10 and the matrices \(M\) and \(N\) are defined as follow: \[\begin{align} M &= (v^\top \otimes I_n)\mathcal{P} = \begin{bmatrix} P^{(1)}v & P^{(2)}v & \dots & P^{(p)}v \end{bmatrix} \in \mathbb{C}^{n\times p},\\ N &= (I \otimes u^\top)\overline{\mathcal{P}} = \begin{bmatrix} (P^{(1)})^* u & (P^{(2)})^* u & \dots & (P^{(p)})^* u \end{bmatrix}\in \mathbb{C}^{n\times p}. \end{align}\]
Proof. We compute the Hessian by differentiating \[\begin{align} d(A+\Delta)v &= \Pi_\mathcal{S}(du\, v^*)v + \Pi_\mathcal{S}(u (dv)^*)v + (A+\Delta)dv,\\ d(A+\Delta)^*u &= \Pi_\mathcal{S}(du\, v^*)^*u + \Pi_\mathcal{S}(u\, (dv)^*)^*u + (A+\Delta)^*\, du. \end{align}\] We recall some identities already used in 16 : for every \(a,b \in \mathbb{C}^n\), we have \[\mathop{\mathrm{vec}}\Pi_\mathcal{S}(ab^*) = \mathcal{P}\mathcal{P^*}(\mathop{\mathrm{vec}}ab^*) = \mathcal{P}\mathcal{P^*} (\overline{b}\otimes a) = \mathcal{P}\mathcal{P^*} (\overline{b}\otimes I_n)a = \mathcal{P}\mathcal{P^*} (I_n\otimes a)\overline{b}.\] These identities lets us simplify a few terms in \(d(A+\Delta)v\) and \(d(A+\Delta)^*u\), using: \[\begin{gather} \Pi_\mathcal{S}(du\, v^*)v = (v^\top \otimes I) \mathop{\mathrm{vec}}\Pi_\mathcal{S}(du\, v^*) = (v^\top \otimes I)\mathcal{P}\mathcal{P^*}(\overline{v}\otimes I_n)du = MM^*du,\\ \Pi_\mathcal{S}(u (dv)^*)v = (v^\top \otimes I) \mathcal{P}\mathcal{P}^*(I_n\otimes u)\overline{dv} = MN^\top \overline{dv},\\ \Pi_\mathcal{S}(u\, (dv)^*)^*u = \mathop{\mathrm{vec}}(u^\top \overline{\Pi_\mathcal{S}(u(dv)^*)}) = (I_n \otimes u^\top)\overline{\mathcal{P}\mathcal{P^*} (I_n\otimes u)\overline{dv}}\\ = (I_n \otimes u^\top)\overline{\mathcal{P}} \mathcal{P}^\top(I_n\otimes \overline{u})dv = NN^*dv,\\ \Pi_\mathcal{S}(du\, v^*)^*u = \mathop{\mathrm{vec}}(u^\top \overline{\Pi_\mathcal{S}(du\, v^*)}) = (I_n\otimes u^\top)\overline{\mathcal{P}\mathcal{P}^*(\overline{v}\otimes I_n) du}\\ = (I_n\otimes u^\top) \overline{\mathcal{P}}\mathcal{P}^\top (v\otimes I_n) \overline{du} = NM^\top \overline{du}. \end{gather}\] ◻
Moreover, for the special case of the projection on a real sparsity structure, this expression can be further simplified.
Lemma 5. Let \(\mathcal{S}\) be the subspace of matrices with sparsity structure \(\mathcal{J}\), i.e., \(B_{ij}=0\) if \((i,j)\not\in\mathcal{J}\). Assume \(A\in\mathcal{S}\) and that \(A,u,v\) have real entries; then, the operator \(H(u,v)\) (over real vectors) is \[\label{Hessian32of32F} H = \begin{bmatrix} K_1 & A+2\Delta\\ (A+2\Delta)^* & K_2 \end{bmatrix},\tag{29}\] where we have set again \(\Delta = \Pi_\mathcal{S}(uv^*)\), and \(K_1, K_2\) are the two diagonal matrices such that \[\begin{align} (K_1)_{ii} = \sum_{j \text{ s.t. }(i,j)\in\mathcal{J}} v_j^2, \quad (K_2)_{jj} = \sum_{i \text{ s.t. }(i,j)\in\mathcal{J}} u_i^2. \end{align}\]
Proof. We may take \((e_ie_j^\top \colon (i,j)\in\mathcal{J})\) as an orthonormal basis of \(\mathcal{S}\). Then, the columns of \(M\) are \(e_ie_j^*v = e_i v_j\), and the columns of \(N\) are \(e_je_i^*u = e_j u_i\). It follows that in Lemma 4 \[\begin{align} MM^* &= MM^\top = \sum_{(i,j)\in\mathcal{J}} e_i v_j v_j e_i^\top = K_1,\\ NN^* &= NN^\top = \sum_{(i,j)\in\mathcal{J}} e_j u_i u_i e_j^\top = K_2,\\ MN^\top &= \sum_{(i,j)\in\mathcal{J}} e_i v_ju_i e_j^\top = \Delta. \end{align}\] ◻
Example 1. Let \(A= \operatorname{diag}(\sigma_1,\sigma_2)\), with \(\sigma_1 > \sigma_2 > 0\), and \(\mathcal{S} = \mathbb{R}^{2\times 2}\) be the trivial sparsity structure satisfied by every \(2\times 2\) matrix (\(\mathcal{J}=\{(1,1),(1,2),(2,1),(2,2)\}\), \(\Pi_\mathcal{S}= I\)). Then, the closest (structured) singular matrix to \(A\) is \(\operatorname{diag}(\sigma_1,0)\), by the Eckhart-Young theorem, corresponding to the rank-1 perturbation \(\Delta = uv^*\) with \(u=- e_2 \sigma_2\), \(v=e_2\). We can verify that \((A+\Delta)v=(A+\Delta)^*u=0\) holds, i.e., 17 is satisfied.
Moreover, in this point \((u,v) = (-\sigma e_2, e_2)\), 29 gives \[H = \nabla^2_{[u,v]}F = \begin{bmatrix} 1 & 0 & \sigma_1 & 0\\ 0 & 1 & 0 & -\sigma_2\\ \sigma_1 & 0 & \sigma_2^2 & 0\\ 0 & -\sigma_2 & 0 & \sigma_2^2 \end{bmatrix}.\] The eigenvalues of \(H\) are the union of those of \[H_1 = \begin{bmatrix} 1 & \sigma_1\\ \sigma_1 & \sigma_2^2 \end{bmatrix}, \quad H_2 = \begin{bmatrix} 1 & -\sigma_2 \\ -\sigma_2 & \sigma_2^2 \end{bmatrix}.\] The matrix \(H_1\) is indefinite, since \(\sigma_1 > \sigma_2\), and hence it has a positive and a negative eigenvalue. The matrix \(H_2\) is semidefinite, with eigenvalues \(0\) and \(1+\sigma_2^2\). The presence of this zero eigenvalue reflects the fact that the solution is overparametrized: \(F(u,v) = F(\alpha u, \frac{1}{\alpha}v)\), and hence \(\frac{d}{d\alpha} F(\alpha u, \frac{1}{\alpha}v) = 0\).
As the Hessian \(H\) is indefinite, the solution of the distance to singularity problem is a saddle point of \(F(u,v)\), not a local minimum or maximum.
4 An algorithm based on the Newton method↩︎
The nonlinear system 17 can be solved with a Newton approach. There are several issues to discuss.
4.1 Fixing the normalization↩︎
Recall that the solutions of 17 are defined up to a normalization factor as in 27 ; and, consequently, the differential \(H\) is singular. In order to avoid a spurious degree of freedom, we wish to compute only solutions with \(\norm{v}=1\). To devise a strategy to enforce this normalization, we first reason on the case in which \(G(u,v)\) is a gradient system. Then, we can choose \(\beta > 0\) and define \[F_{\beta}(u,v) = \frac{1}{2} \norm{A + \Pi_{\mathcal{S}}(uv^*)}_F^2 + \frac{\beta}{4} \left( \norm{v}^2-1 \right)^2.\] The gradient of \(F_{\beta}(u,v)\) is \[G_{\beta}(u,v) = \begin{bmatrix} (A+\Delta)v\\ (A+\Delta)^*u + \beta(\norm{v}^2-1)v \end{bmatrix}, \quad \Delta = \Pi_\mathcal{S}(uv^*).\] Clearly, if \(u,v\) are a stationary point for \(F(u,v)\) and \(\norm{v}=1\) then they are also a stationary point for \(F_{\beta}(u,v)\), since the gradients of both summands are zero.
Even when \(A \not \in \mathcal{S}\) and \(G(u,v)\) is not a gradient system, we can still formulate the modified equation \(G_{\beta}(u,v)=0\): indeed, we can verify directly that any solution \((u,v)\) to \(G(u,v)=0\) in which \(\norm{v}=1\) is also a solution to \(G_{\beta}(u,v)\).
The differential of \(G_{\beta}(u,v)\) is \[H_\beta(u,v) = H(u,v) + \begin{bmatrix} 0 & 0\\ 0 & 2\beta vv^* + \beta(\norm{v}^2-1)I \end{bmatrix},\] where \(H(u,v)\) is defined in Lemma 4. With this differential, we write down the multivariate Newton method for finding a solution \((u,v)\) to the nonlinear system of equations \(G_{\beta}(u,v) = 0\): \[\label{Newton32k43132formula} \begin{bmatrix} u_{k+1}\\ v_{k+1} \end{bmatrix} = \begin{bmatrix} u_{k}\\ v_{k} \end{bmatrix} + \begin{bmatrix} \delta_{u_{k}} \\ \delta_{v_{k}} \end{bmatrix}\tag{30}\] where \[\label{Newton32increment32formula} \begin{bmatrix} \delta_{u_{k}} \\ \delta_{v_{k}} \end{bmatrix} = - H_{\beta}^{-1}(u,v) G_{\beta}(u,v).\tag{31}\]
4.2 Line search↩︎
We modify the update formula 31 to obtain an algorithm with line search, which we formulate as Algorithm 1.
In words, if the Newton step 31 would produce a new iterate \((u_{k+1},v_{k+1}) = (u_k + \delta_{u_k}, v_k + \delta_{v_k})\) that does not reduce the value of \(\norm{G_\beta(u_k,v_k)}\), then we reduce the step-length by half and test a new candidate \(u_{k+1} = u_k + \frac{1}{2}\delta_{u_k}, v_{k+1} = v_k + \frac{1}{2}\delta_{v_k}\); we continue in this fashion halving the step length until the norm is reduced. If \(H_\beta(u,v)\) is positive definite, the Newton direction is always a descent direction, and this strategy should ensure an eventual decrease.
This backtracking strategy is a very simple form of globalization of Newton’s method, to ensure that we produce a decreasing sequence \(\norm{G_\beta(u_k,v_k)}\). We refer the reader to [23] for a broader discussion of line search methods in optimization and more sophisticated strategies.
In practice, we observed that in most cases the choice \(\alpha=1\) (the classical Newton method without line search) is sufficient to get a reduction of the norm, without the need for backtracking steps; but this safeguard increases the robustness of the algorithm with minimal additional cost.
4.3 Starting values↩︎
Since the behavior of the Newton method depends on the choice of the initial point, we describe an effective strategy for selecting initial choices for the vectors \(u,v\) in Algorithm 1.
We recall that in the unstructured case the minimum of 1 is given by \(\Delta = - \sigma_n u_n v_n^*\), where \(u_n,v_n\) are the left and right singular vectors associated with the smallest singular value \(\sigma_n\) of the matrix \(A\). We may use this unstructured minimizer as the starting point of the structured problem, setting \(u_{\text{initial}} = -\sigma_n u_n, \, v_{\text{initial}} = v_n,\) but it turns out that in general there is a better choice than \(\sigma_n\) for the scaling coefficient. In the following, we describe a method to choose \(\sigma\) in an initial value of the form \[\label{eq:initial} u_{\text{initial}} = -\sigma u_n, \quad v_{\text{initial}} = v_n,\tag{32}\] inspired by the ideas in Subsection 2.1.
According to Subsection 2.1, we consider — as a functional to minimize wrt \(E\) (of unit Frobenius norm) — the following: \[\nonumber \widetilde{F}_{\varepsilon}(E) = \sigma_{n} (A + \varepsilon E), \qquad where \qquad \sigma_{n} (A + \varepsilon E) = \sigma_{\min} (A + \varepsilon E).\]
Recall that \(\Delta = \Pi_\mathcal{S}(u v^*) = \varepsilon E\), where \(\norm{E}_F=1\). Let \(\widetilde{G}\) be the gradient of \(E \mapsto \widetilde{F}_{\varepsilon}(E)\) at \(\varepsilon=0\), i.e., \(\widetilde{G}= \Pi_\mathcal{S}(u_n v_n^*)\) with \(u_n\) and \(v_n\) the left and right singular vectors associated to the smallest singular value \(\sigma_{n} (A)\), which we suppose to be simple and nonzero.
We wish to approximate the smallest solution of \[\varphi(\varepsilon):=\widetilde{F}_\varepsilon(E(\varepsilon))=0,\] where \(E(\varepsilon)\) indicates a stationary point of the gradient system 6 , and thus a (local) minimizer of \(\widetilde{F}_\varepsilon(E)\) - for given \(\varepsilon\ge 0\) - over set of matrices of unit Frobenius norm. Clearly at \(\varepsilon=0\), \(E\) does not play any role and \(\varphi(0) = \sigma_n(A)\). In order to compute an effective starting value \(\varepsilon=\varepsilon_0\), we formally apply a Newton step to equation 4 , starting from the value \(\varepsilon=\varepsilon_{-1} = 0\).
For \(k=-1\) with \(\varepsilon_{-1}=0\), the Newton iteration gives \[\label{eq:eps0F} \varepsilon_{\rm initial} := \varepsilon_0 = {0}-\frac{\varphi(0)}{\varphi'(0)} = {}-\frac{\sigma_{n}(A)}{\varphi'(0)}.\tag{33}\] Assuming the smoothness of \(E(\varepsilon)\) wrt \(\varepsilon\), we use [6] to express the derivative of \(\varphi(\varepsilon)\) in explicit form (for \(\varepsilon\) smaller than the structured distance to singularity, i.e. s.t \(\varphi(\varepsilon) > 0\)), namely
\[\varphi'(\varepsilon) = -\| \widetilde{G}(\varepsilon) \|_F = {}-\| \Pi_\mathcal{S}\left(\hat{u}(\varepsilon) v(\varepsilon)^* \right) \|_F,\] where we recall that \(\hat{u} = u/\| u\|\). Therefore at \(\varepsilon=0\), where \(\hat{u}=u_n\) and \(v=v_n\), we have \[\nonumber \varphi'(0) = {}-\| \Pi_\mathcal{S}\left( u_n v_n^* \right) \|_F.\] As a consequence 33 yields \[\label{eq:eps0F2} \varepsilon_{0} = \frac{\sigma_{n}(A)}{\| \Pi_\mathcal{S}\left( u_n v_n^* \right) \|_F};\tag{34}\] we use this value of \(\varepsilon_0\) as the norm for the initial perturbation \(\Delta\). Hence, we choose \(\sigma\) in 32 so that \(\norm{\Delta} = \varepsilon_0\): \[\label{eq:choice-sigma} \sigma= \frac{\varepsilon_0}{\| \Pi_\mathcal{S}\left( u_n v_n^* \right) \|_F} = \frac{\sigma_n}{\| \Pi_\mathcal{S}\left( u_n v_n^* \right) \|^2_F}.\tag{35}\]
Remark 3. The starting value heuristic \(\sigma\) in 35 can be obtained also from a different argument: in an initial value of the form 32 , we select the value of \(\sigma\) that makes \(A+\Delta\) orthogonal (in the Frobenius scalar product) to \(u_n v_n^*\). Indeed, if we plug \(\Delta = \Pi_\mathcal{S}(u_{\text{initial}} v_{\text{initial}}^*) = -\sigma\Pi_\mathcal{S}(u_n v_n^*)\) into \(\langle A+\Delta, u_nv_n^* \rangle=0\) and solve for \(\sigma\), we get \[\sigma = \frac{u_n^*Av_n}{\langle \Pi_\mathcal{S}(u_nv_n^*), u_nv_n^*\rangle} = \frac{\sigma_n}{\norm{\Pi_\mathcal{S}(u_nv_n^*)}_F^2},\] using the fact that \(\Pi_\mathcal{S}\) is an orthogonal projection. Assuming \(A \in \mathcal{S}\), the choice of \(\sigma\) can be also interpreted geometrically as follows: we select \(\sigma\) so that \(A+\Delta\) is orthogonal to the structured steepest descent direction (structured negative gradient) \(\Pi_S(u_n v_n^*)\), i.e., \[\langle A+\Delta,\; \Pi_\mathcal{S}(u_n v_n^*) \rangle = 0.\] Since \(A \in \mathcal{S}\) and \(\Pi_\mathcal{S}\) is the orthogonal projector onto \(\mathcal{S}\), we have \(\langle A,\Pi_\mathcal{S}(u_n v_n^*)\rangle = \langle A,u_n v_n^* \rangle = \sigma_n.\) Therefore, the initial perturbation is obtained by canceling the component of \(A\) along the gradient direction \(\Pi_\mathcal{S}(u_n v_n^*)\).
Multiple initial values↩︎
While the above choice of starting values for \(u\) and \(v\) is reasonable when only a single starting point is used, in some problems it might be necessary to run test several starting values to reduce the risk of getting trapped in a local minimum. Indeed, there may be examples where the nonsingular matrix \(A\) has a set of singular values of small and comparable magnitude. In such settings, adding structure may give an optimization problem 1 with multiple local minima of comparable magnitude; and the initial choice in 32 –35 does not guarantee convergence to the global optimum. Therefore, we suggest running the method in Algorithm 1 for several choices of the initial vectors \(u_{\text{initial}}\), and \(v_{\text{initial}}\). In detail, we select \(K\) pairs of left and right singular vectors \(u_{n-k+1}\) and \(v_{n-k+1}\) associated with the \(n-k+1\)-th smallest singular value \(\sigma_{n-K+1}\), \(k=1,\ldots,K\), and run the method \(K\) times, choosing (for each \(k\)) \[u_{\text{initial}} = \widehat{\sigma}_{n-k+1} u_{n-k+1}, \quad v_{\text{initial}} = v_{n-k+1},\] with \[\label{eeq:choice} \widehat\sigma_{n-k+1} = \frac{\sigma_{n-k+1}}{\norm*{\Pi_\mathcal{S}\left( u_{n-k+1} v_{n-k+1}^* \right)}^2_F}\tag{36}\] The computed final perturbation \(\Delta_* \in \mathcal{S}\), determining the singularity of \(A + \Delta_*\), is then selected as the one of smallest Frobenius norm among the \(K\) computed ones.
An example where the solution benefits from the multiple initializations is presented in the subsequent Section 5.2.
As a cheaper alternative to the one described above, one may select \[u_{\text{initial}} = \widehat{\sigma}_{n-\ell+1} u_{n-\ell+1}, \quad v_{\text{initial}} = v_{n-\ell+1}, \qquad with \qquad \ell = \mathop{\mathrm{arg\,min}}_{1 \le k \le K} \widehat{\sigma}_{n-k+1},\] which would avoid applying the proposed method several times.
5 Numerical experiments↩︎
In this section, we report several numerical experiments that were performed using the Matlab implementation of our method available on https://github.com/fph/NearestSingularAsSystem/. The timings reported refer to a laptop with an Intel Core i5-1135G7 and Matlab R2025b.
5.1 A large-scale matrix↩︎
To illustrate the numerical efficiency of the new method, we take an experiment that appears in [2], [3]: finding the nearest singular matrix to the matrix orani678 in the SuiteSparse Matrix collection, preserving the same sparsity pattern structure. This is an real unsymmetric \(2529\times 2529\) sparse matrix with \(90158\) nonzero elements; it is a case where the minimizer occurs in a rank-drop point for \(M(v)\) in the Riemann-Oracle
method, so regularization is needed. We run Algorithm 1 for this problem, taking advantage of the sparseness: we use svds to compute the left and right singular vectors associated to the smallest singular value, as
a starting point, and we use minres with a very loose tolerance of \(10^{-2}\) to solve the system 31 .
The method takes 5 iterations (with an average of 2083.6 matrix-vector products per iteration inside minres), and converges to a matrix \(A+\Delta\) with \(\sigma_{\min}(A+\Delta)
\approx 1.5481\times 10^{-13}\) and \(\sigma_{\max}(A+\Delta) \approx 3.20\times 10^1\). The algorithm takes less than 3 seconds: this time compares very favorably with the methods in [2] and [3], which take more than 30 seconds to solve the
problem on the same test machine. The computed minima coincide up to at least 7 significant digits.
5.2 A case requiring multiple starting values↩︎
We test the algorithm on the matrix \(C\) available at https://github.com/fph/NearestSingularAsSystem/blob/main/example_starting_value.mat. This is a nonsymmetric \(50\times 50\) matrix with 50% density of nonzeros; it has been obtained as the sparsification of a random-generated orthogonal matrix. This specific matrix has been chosen because the smallest singular value and vectors \(\sigma_n,u_n,v_n\) are not the ones that produce the smallest \(\varepsilon_0\) in 34 . We report in Table 1 the value of the 5 smallest singular values of \(C\), and the magnitude of the local minima \(\norm{\Delta_*}_F\) obtained using them as starting values.
| \(k\) | \(\sigma_k\) | \(\norm{\Delta_*}_F\) |
|---|---|---|
| 46 | 0.1552 | 0.2321 |
| 47 | 0.1015 | 0.1489 |
| 48 | 0.0574 | 0.0793 |
| 49 | 0.0401 | 0.0571 |
| 50 | 0.0389 | 0.0639 |
We see that the smallest norm is obtained with the second-smallest singular value, not the smallest. This example shows that using multiple starting values can produce better solutions.
Since the initial matrix is small, the computation was performed using direct \(O(n^3)\) methods to compute the SVD and to solve the linear system in 31 . The run time of Algorithm 1, repeated 5 times with 5 different starting values, is less than 0.02 seconds.
5.3 An example from polynomial \(\varepsilon\)-GCD computation↩︎
We test the new method also on another example that appears in [2], [20]: given the polynomials \[\label{boito811} p(x) = \gamma_p\prod_{j=1}^{10} (x-\alpha_j), \quad q(x) = \gamma_q\prod_{j=1}^{10} (x-\alpha_j + 10^{-j}), \quad \alpha_j = (-1)^j\tfrac{j}{2},\tag{37}\] we look for the smallest perturbation \(\tilde{p},\tilde{q}\) that gives a pair of polynomials with GCD of prescribed degree \(d\). The normalization coefficients \(\gamma_p\) and \(\gamma_q\) appearing in 37 are chosen so that the vectors of coefficients of \(p\) and \(q\) have Euclidean norm 1.
This problem can be formulated as the distance to singularity problem for a Sylvester matrix \(A_S\) defined as \[A_S = \begin{bmatrix} \frac{1}{\sqrt{\deg(q)-d+1}} \mathcal{T}_{p} & \frac{1}{\sqrt{\deg(p)-d+1}} \mathcal{T}_{q} \end{bmatrix},\] composed of two blocks with Toeplitz structure \[\mathcal{T}_{p} = \underbrace{\begin{bmatrix} p_0\\ p_1 & p_0\\ \vdots & p_1 & \ddots\\ p_{10} & \vdots & \ddots & p_0\\ & p_{10} & \ddots & p_1\\ & & \ddots & \vdots\\ & & & p_{10}\\ \end{bmatrix}}_{\text{\deg(q)-d+1 columns}}, \quad \mathcal{T}_{q} = \underbrace{\begin{bmatrix} q_0\\ q_1 & q_0\\ \vdots & q_1 & \ddots\\ q_{10} & \vdots & \ddots & q_0\\ & q_{10} & \ddots & q_1\\ & & \ddots & \vdots\\ & & & q_{10}\\ \end{bmatrix}}_{\text{\deg(p)-d+1 columns}}.\] We refer to [2] for more details. The problem is quite challenging numerically, as the objective function \(\norm{\Delta}_F^2\) is smaller than the machine precision \(\approx 10^{-16}\) for \(d \leq 5\).
We report the results obtained by the new method, and a comparison with [2], [20], in Table 2.
| Approximate GCDs of [boito811] | |||||
| \(d\) | Algorithm [algo:main] | Riemann-Oracle | \(VP_S\) [20] | \(VP_g\), \(VP_h\) [20] | UVGCD |
| 9 | \(3.9964\cdot 10^{-3}\) | \(3.9964\cdot 10^{-3}\) | \(4.00\cdot10^{-3}\) | \(4.00\cdot10^{-3}\) | \(3.996\cdot 10^{-3}\) |
| 8 | \(1.7288\cdot 10^{-4}\) | \(1.7288\cdot 10^{-4}\) | \(1.73\cdot 10^{-4}\) | \(1.73\cdot 10^{-4}\) | \(1.729\cdot 10^{-4}\) |
| 7 | \(7.0890\cdot 10^{-6}\) | \(7.0890\cdot 10^{-6}\) | \(\underline{3.93\cdot 10^{-4}}\) | \(7.09\cdot 10^{-6}\) | \(7.089\cdot 10^{-6}\) |
| 6 | \(1.8293\cdot 10^{-7}\) | \(1.8293\cdot 10^{-7}\) | \(\underline{2.7\cdot 10^{-5}}\) | \(1.83\cdot 10^{-7}\) | \(1.829\cdot 10^{-7}\) |
| 5 | \(\underline{2.0201}\cdot 10^{-9}\)(*) | \(4.4\underline{913}\cdot 10^{-9}\) | \(\underline{1.72\cdot 10^{-5}}\) | \(4.49\cdot 10^{-9}\) | \(4.487\cdot 10^{-9}\) |
| 4 | \(\underline{2.7776}\cdot 10^{-11}\) (*) | \(1.8293\cdot 10^{-7}\) | \(\underline{1.55\cdot 10^{-14}}\) | \(8.40\cdot 10^{-11}\) | \(8.40\cdot 10^{-11}\) |
The results show that the method obtains reliable results up to the point where \(\norm{\Delta}_F^2\) becomes smaller than the machine precision. The comparison with the results in [20] is interesting: the new method performs better than SLRA on the same matrix formulation \(A_S\) of the polynomial GCD problem (column \(VP_S\) in [20]), even though SLRA obtains better results when run on a different matrix formulation of the problem (columns \(VP_g\), \(VP_h\)) that does not fall in the exact same framework described in Section 2.2. The results of \(VP_g\), \(VP_h\) from [20] match those of a state-of-the-art specialized algorithm (column \(UVGCD\), reported with more significant digits from [24]). Overall, the results suggest that the new solution algorithm is reliable up to the point where \(\norm{\Delta}_F^2\) becomes smaller than the machine precision: it is as accurate as the one in [2] and slightly more than the one in [20], when run on the same matrix \(A_S\).
6 Conclusions↩︎
We have introduced and studied a new method to numerically approximate the structured distance to singularity. After a theoretical comparison between two existing methods, we propose a novel approach, which is independent of the existing ones. The technique solves the nearness problem via a Newton method on a system of nonlinear equations, avoiding the nested optimizations arising in [2] and [3]. We applied this technique both to large, sparse matrices and smaller problems with Toeplitz-like structures. Further research directions include expanding the theoretical analysis of the proposed method, and generalizing it to distance to stability problems.
6.0.0.1 Acknowledgements
All the authors are affiliated to the Italian INdAM-GNCS (Gruppo Nazionale di Calcolo Scientifico). FP acknowledges the support by the National Centre for HPC, Big Data and Quantum Computing–HPC, CUP B83C22002940006, funded by the European Union–NextGenerationEU, and by the Italian Ministry of University and Research through the PRIN project 2022 “MOLE: Manifold constrained Optimization and LEarning”, CUP B53C24006410006. SS acknowledges the support by the European Union (ERC consolidator, eLinoR, no 101085607). NG acknowledges that his research was supported by funds from the Italian MUR (Ministero dell’Università e della Ricerca) within the PRIN 2022 Project “Advanced numerical methods for time dependent parametric partial differential equations with applications” and the PRIN-PNRR Project “FIN4GEO”. He also acknowledges funding from the Dipartimento di Eccellenza 2023–2027 project awarded to the Gran Sasso Science Institute (GSSI) by the Italian Ministry of University and Research (MUR).
References↩︎
Istituto di Scienza e Tecnologie dell’Informazione “Alessandro Faedo”, CNR, Pisa, Italy.
email: miryam.gnazzo@isti.cnr.it↩︎Division of Mathematics, Gran Sasso Science Institute, L’Aquila, Italy.
email: nicola.guglielmi@gssi.it↩︎Department of Computer Science, University of Pisa, Pisa , Italy.
email: federico.poloni@unipi.it↩︎Department of Mathematics and Operational Research, University of Mons, Mons, Belgium.
email: Stefano.SICILIA@umons.ac.be↩︎