Manifold–affine intersection problems, structured low-rank approximation, regularization, inexact projection, intrinsic transversality, transversality, local convergence

1 Introduction↩︎

Structured low-rank approximation (SLRA) seeks to approximate a given structured matrix (\(M\)) by a matrix (\(M_r\)) with the same structure and prescribed rank. We use SLRA mainly as a motivating setting for local manifold–affine intersection problems. Typical structures include Hankel, Toeplitz, Sylvester, and more general affine matrix structures [1], [2]. Unlike classical low-rank approximation, which is solved explicitly by the SVD, SLRA must enforce affine structure and rank constraints simultaneously. It arises in signal processing, system identification, and computer algebra [1], [3], [4]; low-rank Hankel approximation is related to exponential data fitting, frequency estimation, and denoising, while Sylvester-structured low-rank approximation is important in approximate polynomial GCD computations [1], [4], [5]. Geometrically, SLRA is a representative manifold–affine intersection problem, with the structure encoded by an affine subspace and the rank constraint by a smooth fixed-rank manifold.

Let \(\mathcal{M}_{p,q}(\mathbb{R})\) denote the space of real \(p\times q\) matrices, equipped with the Frobenius inner product \(\langle M_1,M_2\rangle=\operatorname{tr}(M_1^\top M_2)\) and norm \(\|M\|=\sqrt{\operatorname{tr}(M^\top M)}\). Denote by \(\mathcal{D}_r\) the subset of \(\mathcal{M}_{p,q}(\mathbb{R})\) consisting of matrices of rank \(r\). In the SLRA setting, all local statements involving \(\mathcal{D}_r\) are understood in a neighborhood of a rank-\(r\) point \(x^*\in E\cap\mathcal{D}_r\) with \(\sigma_r(x^*)>0\), where \(\mathcal{D}_r\) is a smooth fixed-rank manifold and the metric projection onto \(\mathcal{D}_r\) is locally single-valued. The local SLRA problem considered in this paper can be formulated as follows.

Computing a metric projection onto \(E\cap \mathcal{D}_r\) is generally difficult. In practice, one usually seeks a nearby feasible point in \(E\cap \mathcal{D}_r\), which can be regarded as a local approximation to the projection onto the manifold–affine intersection. Many algorithms have been proposed for this purpose. Among them, the Cadzow algorithm [6], which alternates between the affine structure set and the fixed-rank manifold, is simple and robust, and converges locally linearly under separability-type intersection conditions [7], [8]. To obtain faster local convergence, Newton-SLRA [9] replaces the fixed-rank manifold by its tangent space at each iteration and achieves quadratic convergence under transversality. Other approaches include structured total least squares methods [10], [11], Riemannian optimization [12], [13], alternating minimization [14], variable projection [15], gradient systems [16], and algebraic methods [17].

The quadratic convergence of Newton-SLRA, however, relies on a restrictive geometric requirement: in the SLRA setting, the affine subspace \(E\) and the fixed-rank manifold \(\mathcal{D}_r\) must intersect transversally. When transversality fails, the tangent-space intersection in the Newton-SLRA update may be empty, and the associated linear system may become singular or severely ill-conditioned. Thus, although alternating projections remain applicable under weaker intersection conditions, the fast local behavior of Newton-SLRA is no longer available. This motivates the construction of a Newton-type method that is well defined under weaker geometry and recovers fast convergence under transversality.

To this end, we propose Regularized Newton-SLRA (RN-SLRA) for local manifold–affine intersection problems, with SLRA as a representative application. RN-SLRA replaces the exact tangent-space intersection step in Newton-SLRA by a regularized subproblem. Under intrinsic transversality, we prove local linear convergence of RN-SLRA and also record a local equivalence between intrinsic transversality and the two directional versions of separability in the \(C^2\) manifold–affine setting. Under transversality, we obtain faster local convergence by choosing \[\mu_k=c\,r_k^\rho,\qquad r_k:=d_{M}(x_k)=\|x_k-P_{M}(x_k)\|,\qquad c>0,\quad \rho\in[0,1].\] Since \(r_k\to0\), the regularization increasingly enforces the tangent-space condition underlying Newton-SLRA. With this choice, RN-SLRA converges locally with order at least \(1+\rho\), and in particular quadratically when \(\rho=1\). We call the method geometry-adaptive because the regularization parameter is chosen according to the manifold residual \(d_M(x_k)\), so that the step interpolates between alternating-projection behavior away from the intersection and Newton-SLRA behavior near the intersection.

For large-scale SLRA problems, the projection onto the fixed-rank manifold may itself be computationally expensive. We therefore further introduce Inexact Regularized Newton-SLRA (iRN-SLRA), which replaces the exact manifold projection by a \(\sigma\)-quasioptimal projection [18]. The inexact variant preserves the local convergence rates of RN-SLRA under the corresponding intersection assumptions, while allowing cheaper projection computations. The local convergence properties of the methods discussed above are summarized in Table 1.

We evaluate RN-SLRA and iRN-SLRA on constructed degenerate SLRA instances and low-rank Hankel approximation problems at both small and large scales. The numerical results are consistent with the theoretical motivation above: on instances where transversality fails, RN-SLRA remains convergent, whereas Newton-SLRA may diverge. In addition, iRN-SLRA exhibits convergence behavior comparable to that of RN-SLRA while reducing the computational cost of the manifold projection step, with more pronounced savings on large-scale problems.

2 Preliminaries↩︎

We consider the general problem of approximating the projection of a point \(x\) onto the intersection \(X:=M\cap E\) of a smooth embedded manifold \(M\) and an affine subspace \(E\) in a Euclidean space, assuming without loss of generality that \(M\neq E\); in particular, the SLRA problem is recovered when \(M=\mathcal{D}_r\).

2.1 Notations and Basic Facts↩︎

We use \(d_C(x)\) to denote the distance from \(x\) to the set \(C\), and \(P_C(x)\) to represent the projection of \(x\) onto the set \(C\). \(B(x,r)\) denotes the ball centered at \(x\) with radius \(r\). We use \(\hat{u} = \frac{u}{\|u\|}\) to represent the unit vector of \(u\). For any affine subspace \(A\), let \(A^0\) denote its underlying vector space, so that \(A=x+A^0 \text{ for any }x\in A.\) In particular, we write \(E=x+L \text{ for any }x\in E,\) where \(L:=E^0\). For an affine manifold \(V\) and a point \(x\in V\), we denote by \(T_xV\) the affine tangent space of \(V\) at \(x\), and by \(N_xV\) the affine normal space of \(V\) at \(x\). Their underlying vector spaces are denoted by \(T_xV^0\) and \(N_xV^0\). \(L^{\perp}\) denotes the orthogonal complement of \(L\).

We first recall some intersection conditions, including the transversality condition [7], [19], [20], the clean intersection condition [21], the intrinsic transversality condition [22], [23], the separability condition [24], and the subtransversality condition [25].

We have the following lemma under transversality:

This lemma follows from [9]. A proof can be found in [26].

The notion of separability was first introduced in [24] under the term 0-separability. This property is particularly suitable for the analysis of alternating projections and inexact alternating projection methods [8].

In our work, intrinsic transversality serves as the main geometric assumption for the local linear convergence analysis.

Under Assumption 1, intrinsic transversality is equivalent to separability in either order; see Lemma 2 below. Although separability is not symmetric in general, in the present \(C^2\) manifold–affine setting the two directional versions are equivalent. More precisely, \(E\) intersects \(M\) separably at \(\bar x\) if and only if \(M\) intersects \(E\) separably at \(\bar x\), and these are both equivalent to intrinsic transversality.

Under the present \(C^2\) embedded manifold–affine setting and with the definitions used here, the following relationships hold locally, rather than as unconditional statements for arbitrary sets: transversality implies clean intersection; clean intersection is equivalent to intrinsic transversality[27]; intrinsic transversality implies subtransversality; and intrinsic transversality is equivalent to separability in either order. See [23], [24], [27].

2.2 Cadzow Algorithm↩︎

The alternating projection (AP) method has a long history [28], [29]. Recent studies have analyzed AP and inexact AP methods for various sets [30]–[32], together with their local linear convergence properties [8], [24], [31], [33]. In the context of structured low-rank approximation, the Cadzow algorithm [6], [34] can be viewed as an alternating projection method between the affine space \(E\) and the low-rank manifold \(M\). More precisely, starting from \(x_k\in E\), one step of the Cadzow iteration is given by \[\begin{align} x_{k+1}=P_E P_M(x_k). \end{align}\]

In the local rank-\(r\) neighborhood considered in this paper, the projection onto \(\mathcal{D}_r\) is computed by the truncated SVD, as a local consequence of the Eckart–Young theorem:

Then Cadzow algorithm can be formulated as: \[\begin{align} &U,S,V \leftarrow \text{SVD}(x_k),\qquad y_k=U_rS_rV_r^\top, \\ &x_{k+1}=x_k+\sum_{i=1}^d \langle y_k-x_k,E_i \rangle E_i. \end{align}\] where \(U_r\) denotes the matrix consisting of the first \(r\) columns of \(U\), \(V_r\) denotes the matrix consisting of the first \(r\) columns of \(V\), \(S_r\) denotes the \(r\times r\) top-left submatrix of \(S\), \(d\) is the dimension of \(E\), and \(\{E_i\}_{i=1}^d\) denotes an orthonormal basis of \(L\). The update formula for \(x_{k+1}\) is obtained by \[x_{k+1}=P_E(y_k)=x_k+P_L(y_k-x_k).\]

If \(E\) intersects \(M\) transversally at \(\bar x\in X\) and the initial point is sufficiently close to \(\bar x\), the Cadzow algorithm will converge locally linearly to a point \(x_{\infty}\in X\) [7]. In fact, local linear convergence can already be guaranteed when \(E\) intersects \(M\) separably [8].

2.3 Newton-SLRA Algorithm↩︎

The Newton-SLRA algorithm is motivated by Newton’s method. In each iteration of the classical Newton-SLRA, given the current iterate \(x_k \in E\), one first computes its projection onto the manifold \(y_k = P_{M}(x_k)\), and then solves a linear least-squares problem: \(x_{k+1} = P_{E \cap T_{y_k}M}(x_k),\) i.e., the closest point to \(x_k\) in the intersection of \(E\) and the tangent space \(T_{y_k}M\). Under the transversality condition, this algorithm converges locally quadratically. When \(M=\mathcal{D}_r\), the algorithm can be formulated as: \[\label{eq:newton-slra} \begin{align} &U,S,V\leftarrow \mathrm{SVD}(x_k),\qquad y_k = U_r S_r V_r^\top,\\ &N_{ij}\leftarrow u_i v_j^\top,\quad i=r+1,\dots,p,\;\;j=r+1,\dots,q,\\ &A\leftarrow \big(\langle N_{ij},E_\ell\rangle\big)_{(i,j),\ell}, \qquad b\leftarrow \big(\langle N_{ij},\,y_k-x_k\rangle\big)_{(i,j)},\\ &x_{k+1}=x_k+\displaystyle\sum_{\ell=1}^d (A^\dagger b)_\ell E_\ell. \end{align}\tag{1}\]

However, when transversality fails, the matrix \(A\) in 1 may lose rank at the solution. In the fixed-rank SLRA setting, full row rank of the tangent-space linear system is equivalent to transversality of \(E\) and \(\mathcal{D}_r\) at the intersection point. Hence, near a nontransversal solution, the associated least-squares problem may become singular or severely ill-conditioned.

Another issue is that the well-definedness of Newton-SLRA relies on the non-emptiness of \(E \cap T_{y_k} M\). Under transversality, there exists a positive lower bound \(\alpha\) for the angle between \(E\) and \(T_{y_k} M\) in some neighborhood of \(\bar{x}\), but this cannot be guaranteed under intrinsic transversality, as \(E \cap T_{y_k} M\) may be empty.

These issues inspired us to modify the Newton-SLRA algorithm to design a new algorithm, called Regularized Newton-SLRA (RN-SLRA), which is able to handle the above issues.

3 Regularized Newton-SLRA↩︎

We introduce a regularization term and define the next iterate as: \[\label{eq:rnslra} \begin{align} &y_k=P_M(x_k),\\ &x_{k+1} = \arg\min_{x \in E} \left\{ \dfrac{1}{2} \| x - y_k \|^2 + \dfrac{1}{2\mu_k} \| P_{N_{y_k}M}(x - y_k) \|^2 \right\}, \quad \mu_k > 0. \end{align}\tag{2}\] where \(P_{N_{y_k}M}\) is the orthogonal projector onto the normal space \(N_{y_k}M\), and \(\mu_k\) is a regularization parameter. In the algorithm studied below, \(\mu_k\) is chosen adaptively according to: \[\mu_k = c \cdot r(x_k)^\rho,\quad r(x_k)=\|x_k-y_k\|=d_M(x_k).\] where \(c>0\) is a constant, \(0\le\rho\le 1\) is a parameter. This choice of the regularization parameter is inspired by the residual-dependent regularization strategy used in the Inexact Regularized Proximal Newton method (IRPN) [35], and is also related to earlier regularized Newton methods for convex minimization with singular solutions [36] as well as improved proximal Newton implementations such as newGLMNET [37].

Since \(E\) is affine, the objective function in 2 for solving \(x_{k+1}\) is a strongly convex quadratic problem with a unique solution, which can be obtained by solving a linear system. Thus, the well-definedness of the RN-SLRA step does not require the transversality condition between \(E\) and \(M\).

3.1 Derivation of the Linear System↩︎

Let \(\{E_1,\dots,E_d\}\) be an orthonormal basis of \(L\), set \(b=x_k-y_k\), and write \(N=N_{y_k}M\). Choose an orthonormal basis \(\{N_j\}_{j=1}^s\) of \(N\). Since \(E=x_k+L\), any \(x\in E\) can be written as \(x=y_k+b+\sum_{i=1}^d\alpha_iE_i\). Substituting this expression into 2 and using the first-order optimality condition gives \[\begin{align} \label{eq:sys} \left(I_d+\frac{1}{\mu_k}A^\top A\right)\alpha = -\beta-\frac{1}{\mu_k}A^\top \eta, \end{align}\tag{3}\] where \(A_{ji}=\langle N_j,E_i\rangle, \eta_j=\langle N_j,b\rangle,\) and \(\beta_i=\langle b,E_i\rangle.\) Since \(A^\top A\succeq0\) and \(\mu_k>0\), the coefficient matrix is positive definite. The next iterate is \(x_{k+1}=x_k+\sum_{i=1}^d\alpha_iE_i\).

For the fixed-rank manifold \(M=\mathcal{D}_r\), this yields Algorithm 1.

3.2 Limiting Interpretation of the Regularized Step↩︎

Assume that \(E\) and \(M\) intersect transversally at \(\bar x\), so that \(E\cap T_{y_k}M\) is nonempty for \(y_k\) sufficiently close to \(\bar x\). The role of \(\mu_k\) can be understood from two limiting regimes.

If \(\mu_k\to\infty\), the penalty term vanishes and the subproblem reduces to \(\min_{x\in E}\|x-y_k\|^2\). Hence the step becomes \(p_k=P_E(y_k)\), namely the alternating projection step. If \(\mu_k \to 0\), the penalty term enforces \(P_{N_{y_k}M}(x-y_k)=0\), or equivalently \(x\in T_{y_k}M\). Thus the limiting problem is \(\min_{x\in E\cap T_{y_k}M} \|x-y_k\|^2 .\) Since \(y_k=P_M(x_k)\), we have \(x_k-y_k\in N^0_{y_k}M\). Hence, for every \(x\in T_{y_k}M\), \(\|x-x_k\|^2 = \|x-y_k\|^2+\|x_k-y_k\|^2 .\) Therefore, minimizing the distance to \(y_k\) over \(E\cap T_{y_k}M\) is equivalent to minimizing the distance to \(x_k\) over the same set. The limiting step is the Newton–SLRA point \(q_k=P_{E\cap T_{y_k}M}(x_k),\) provided that \(E\cap T_{y_k}M\) is nonempty.

Thus, for \(\mu_k\in(0,\infty)\), the regularized step interpolates between alternating projections and Newton-SLRA: smaller \(\mu_k\) gives a more Newton-like step, while larger \(\mu_k\) gives a more projection-like step. This interpretation relies on the local nonemptiness of \(E\cap T_{y_k}M\), which follows from transversality but may fail under Assumption 2; it is used only as intuition and not in the convergence analysis below.

3.3 Inexact Regularized Newton-SLRA↩︎

We replace \(y_k=P_M(x_k)\) in 2 by a quasioptimal projection and obtain iRN-SLRA: \[\label{eq:iRN-SLRA} \begin{align} &y_k\in M \;\text{such that}\;\|x_k-y_k\|\le \sigma\, d_M(x_k), \quad \sigma\ge 1, \\ &x_{k+1}=\displaystyle \arg\min_{x\in E} \left\{ \frac{1}{2}\|x-y_k\|^2 +\frac{1}{2\mu_k}\bigl\|P_{N_{y_k}M}(x-y_k)\bigr\|^2 \right\}, \quad \mu_k>0. \end{align}\tag{4}\] When \(\sigma=1\), iRN-SLRA reduces to RN-SLRA. For \(\sigma>1\), \(y_k\) is a \(\sigma\)-quasioptimal projection of \(x_k\) onto \(M\), in the sense of [18]. This quasioptimal model is well suited to low-rank approximation: truncated SVD gives the exact projection with \(\sigma=1\), while randomized SVD provides a practical realization in expectation or with high probability under standard randomized low-rank approximation bounds [38]; partial SVD can also be used in practice [39]. Compared with a full SVD, these approaches can significantly reduce the computational cost, especially for large-scale problems.

For an \(m\times n\) rank-\(r\) SLRA problem with affine dimension \(d\), one RN-SLRA iteration consists of: (i) computing the rank-\(r\) projection, implemented by a truncated or full SVD; (ii) forming the normal-coordinate matrix \(A_k\in\mathbb{R}^{(m-r)(n-r)\times d}\); and (iii) solving a \(d\times d\) symmetric positive definite linear system. Thus the dominant cost is typically the rank-\(r\) projection and the construction of \(A_k\), while the \(d\times d\) solve is moderate when \(d\) is small or structured. For Hankel/Toeplitz structures, the affine basis and matrix-vector products can often be represented implicitly, and randomized or partial SVD can reduce the projection cost in large-scale settings.

4 Local Convergence Analysis↩︎

In this section, we establish the local convergence properties of the proposed algorithms. We first prove that RN-SLRA converges linearly under intrinsic transversality. We then show that, under the stronger transversality condition, the convergence rate improves to quadratic convergence. Finally, we extend the analysis to the inexact variant iRN-SLRA. Throughout this section, \(\{x_k\}\) denotes the sequence of iterates generated by the corresponding algorithm.

4.1 Linear Convergence of RN-SLRA under Intrinsic Transversality↩︎

Both of the above lemmas can be derived from subtransversality, which itself follows from intrinsic transversality. These two lemmas show that there exists a neighborhood \(U\) of \(\bar x\) such that the distance from a point in \(E\cap U\) to \(M\), or from a point in \(M\cap U\) to \(E\), is of the same order as its distance to \(X\).

In the above proof, we introduced the operator \(B_k=I+\frac{1}{\mu_k}P_{N_k}P_L: N_k\to N_k\) in 16 , and established the bound \(\|B_k^{-1}\|\le 1\). We emphasize that a sharper estimate requires a uniform lower bound on the action of \(P_L\) on \(N_k\). More precisely, under transversality one has, after shrinking the neighborhood if necessary, there exists \(\eta>0\) such that \(\|P_L n\|\ge \eta\|n\|\) for all \(n\in N_k\). Equivalently, \(P_{N_k}P_L|_{N_k}\) is uniformly positive definite on \(N_k\). We will use this stronger estimate in Theorem 3.

The following result is a direct corollary of the previous lemma.

The next result shows that the limit point selected by Algorithm RN-SLRA has the same order of accuracy as the metric projection onto \(X\). After shrinking the neighborhood if necessary, \(X\) is a \(C^2\) embedded submanifold and the local projection \(P_X\) is single-valued.

Before the next remark, we recall that intrinsic transversality is equivalent to clean intersection for intersections of \(C^2\) embedded submanifolds [27]. Therefore, in the present manifold–affine setting, Theorem 2 can be interpreted within the clean intersection framework.

4.2 High-order Convergence of RN-SLRA under Transversality↩︎

We have proved that RN-SLRA converges linearly to a point in \(X\) under the assumption that \(E\) intersects \(M\) intrinsically transversally at \(\bar x\in X\). If we strengthen the intersection condition to transversality, we obtain sharper local estimates for RN-SLRA. These estimates give a linear-type bound when \(\rho=0\), superlinear convergence when \(0<\rho<1\), and quadratic convergence when \(\rho=1\). Since transversality implies intrinsic transversality, the lemmas we have proved before still hold.

4.3 Convergence Analysis of iRN-SLRA↩︎

In this subsection, we show that iRN-SLRA still enjoys local linear convergence under the intrinsic transversality condition. Under the transversality condition, it also retains the local superlinear and quadratic convergence properties.

We use the following notation throughout. Let \[\tilde{r}_k:=\|x_k-y_k\|, \qquad \mu_k:=c\,\tilde{r}_k^\rho, \quad c>0,\;\rho\in[0,1].\] We also keep the notation \(r_k:=d_M(x_k)\). Furthermore, we denote by \(p_k:=P_E(y_k)\) the projection of \(y_k\) onto \(E\).

We begin by relating the exact residual \(r_k\) and the inexact residual \(\tilde{r}_k\).

The next estimate is the inexact analogue of the step-size bound proved earlier for RN-SLRA. Its proof is identical, once \(r_k\) is replaced by \(\tilde{r}_k\).

We now prove local linear convergence under intrinsic transversality. The argument reuses the residual-contraction proof for RN-SLRA, with the alternating-projection comparison point \(p_k=P_E(y_k)\).

Under transversality, the proof simplifies further: the comparison point \(p_k\) still enters only through first-order information, while the regularized Newton correction suppresses the normal component by an additional factor \(\mu_k\). This yields the same high-order estimate as in the exact-projection case.

5 Experimental Results↩︎

In this section, we test our algorithm on some applications of structured low-rank approximation, and compare it with representative existing methods. ² To demonstrate the ability of our algorithm to handle problems with poor intersection conditions, we will construct some specific examples. All experiments were conducted in Python 3.9.13 on an Ubuntu 22.04 server equipped with dual AMD EPYC 9754 128-core processors (512 logical threads) and 1 TB of RAM.

5.1 Small-scale SLRA Problems↩︎

In this section, we construct small-scale SLRA problems that do not satisfy the transversality condition, in order to illustrate the following phenomena:

• RN-SLRA may exhibit either superlinear or only linear local convergence, depending on the geometry;

• In some cases, Newton-SLRA may diverge while RN-SLRA converges;

• Under intrinsic transversality, the theory guarantees only local linear convergence, but RN-SLRA can still converge much faster than Cadzow in ill-conditioned examples.

All small-scale examples are SLRA problems with artificially constructed affine subspaces and the \(4\times4\) rank-\(2\) manifold \(\mathcal{D}_2\). They are constructed near \(x^*=\mathop{\mathrm{diag}}(4,\sigma_2,0,0)\), where \(\sigma_2>0\) is small. The parameter \(\sigma_2\) controls the local geometry of the rank constraint: smaller \(\sigma_2\) leads to a more ill-conditioned low-rank geometry and worsens the conditioning of the tangent-space system [44], [45].

We choose the affine subspace in the form \(E=x^*+\operatorname{span}(B_1,B_2,B_3,B_4),\) where the directions \(B_i\) are obtained by orthonormalizing prescribed raw directions. These raw directions are chosen mainly in the tangent space \(T_{x^*}\mathcal{D}_r\), with a small number of normal components added. Thus each direction can be decomposed as \(B_i=T_i+N_i,\) where \(T_i\in T_{x^*}\mathcal{D}_r,\) and \(N_i\in N_{x^*}\mathcal{D}_r .\)

For the diagonal point \(x^*=\mathop{\mathrm{diag}}(4,\sigma_2,0,0)\), we use the natural singular-vector coordinates. In these coordinates, the entries coupling the leading and trailing two coordinates give tangent directions, whereas the lower-right \(2\times2\) block gives normal directions. For instance, \(E_{13},E_{31},E_{24},E_{42}\in T_{x^*}\mathcal{D}_r\), while \(E_{33},E_{34},E_{43},E_{44}\in N_{x^*}\mathcal{D}_r\).

The normal components are chosen to be linearly dependent. Indeed, by the definition of the matrix \(A\) in 3 , for a normal basis element \(N_{ab}=u_av_b^\top\), we have \(A_{(a,b),i}(x^*) =\langle N_{ab},B_i\rangle =u_a^\top B_i v_b =u_a^\top N_i v_b .\) Hence \(A(x^*)\) depends only on the normal components \(N_i\). If these components are linearly dependent, then \(A(x^*)\) is rank deficient. Since transversality at \(x^*\) is equivalent to the full-rank condition of \(A(x^*)\), the constructed examples violate the transversality condition.

Finally, the initial point is taken as \(x_0=x^*+\sum_{i=1}^4 a_i B_i,\) which is a small perturbation of \(x^*\) within the affine space \(E\). The explicit choices of the raw directions and coefficients \(a_i\) are given in the individual examples below.

5.1.1 Example 1: Newton-SLRA diverges while RN-SLRA converges↩︎

We first take \(x^* = \operatorname{diag}(4,0.02,0,0),\) and define directions: \[\begin{cases} \widetilde{B}_1 = E_{13}+0.7E_{31}+E_{33},& \widetilde{B}_2 = E_{24}-0.6E_{42}+E_{33},\\ \widetilde{B}_3 = E_{34}+0.2E_{12},& \widetilde{B}_4 = E_{44}-0.1E_{21}.\\ \end{cases}\]

The actual basis \(B_1,\dots,B_4\) is obtained from \(\widetilde{B}_1,\dots,\widetilde{B}_4\) by Gram-Schmidt orthonormalization. The initial point is then fixed as \(x_0 = x^* + 0.02\sum_{i=1}^4 a_i B_i,\) with \(a = (-0.98912135,\,-0.36778665,\,1.28792526,\,0.19397442).\)

We compare the numerical performance of Newton-SLRA and Algorithm 1 (denoted as RN-SLRA). We set the parameters of RN-SLRA to \(c=3\times 10^{-7}\) and \(\rho=1\). Our numerical results are presented in Fig. 2.

RN-SLRA converges in \(7\) steps, taking \(0.001852s\), while Newton-SLRA diverges. From the right plot, we can see that the condition number of matrix \(A\) becomes extremely large during the iteration process of RN-SLRA, while RN-SLRA still converges rapidly. This example illustrates the stabilizing effect of the regularization term in a nontransversal SLRA problem.

To further check that the above behavior is not caused by a single specially chosen initial point, we also perform a random initialization test on the same example. We generate \(x_0 = x^* + \alpha \sum_{i=1}^4 a_i B_i,\) where \(a=(a_1,\ldots,a_4)\) is obtained by sampling a standard Gaussian vector in \(\mathbb{R}^4\) and normalizing it to satisfy \(\|a\|_2=1\). The radius \(\alpha\) is sampled independently from the uniform distribution \(\alpha \sim U(0.01,0.1).\) We use \(1000\) random initial points. The parameters of RN-SLRA are kept the same as above. The success rates are reported in Table 2.

The result shows that Newton-SLRA is sensitive to the initialization in this nontransversal example, whereas RN-SLRA succeeds for all tested initial points. This supports the observation from Fig. 2 that the regularization stabilizes the Newton iteration near this degenerate intersection.

5.1.2 Example 2: linear convergence of RN-SLRA and comparison with Cadzow↩︎

In the previous example, we observed the superlinear convergence of RN-SLRA. We now construct an example for which RN-SLRA exhibits only linear local convergence. This shows that the theoretical linear convergence rate of RN-SLRA is, in general, not improvable. Moreover, using the same example, we compare RN-SLRA with the Cadzow algorithm and show that RN-SLRA can converge much faster, although both algorithms are linearly convergent.

Similar to the previous construction process, we take \(x^*=\operatorname{diag}(4,0.1,0,0),\) and \[\begin{cases} \widetilde{B}_1 = E_{13}+0.6E_{31}+1.0E_{33},& \widetilde{B}_2 = E_{24}-0.5E_{42}+1.0E_{33},\\ \widetilde{B}_3 = E_{34}+0.1E_{12}+0.2E_{44},& \widetilde{B}_4 = E_{43}-0.1E_{21}+0.6E_{44}+0.4E_{33}. \end{cases}\] The initial point is chosen as \(x_0=x^*+0.01\sum_{i=1}^4 a_i B_i,\) with \(a=(1,-1,0.5,0.2).\) We set the parameters of RN-SLRA to \(c=0.001\) and \(\rho=1\). On this example, RN-SLRA converges linearly in \(10\) steps, taking \(0.002545\)s, with an approximate convergence rate \(q=0.2501\).

We further compare the numerical performance of RN-SLRA and Cadzow on the same example with the same parameters. For both algorithms, the maximal number of iterations is set to \(1000\). The numerical results are presented in Table 3. The residual ratio is the median of the last five values of \(r_{k+1}/r_k\), with \(r_k=\|x_k-P_{D_r}(x_k)\|\).

The results show that RN-SLRA converges much faster than Cadzow on this example. This behavior is consistent with the local geometry of the problem: alternating projections can be slow when the two sets meet with a very small effective angle, as in the simple case of two nearly parallel intersecting lines. Our constructed example exhibits a similar local geometry near the intersection point. Thus, this example illustrates that RN-SLRA may have only linear local convergence under nontransversality, but can still be much faster than Cadzow.

5.2 Large-scale SLRA Problems↩︎

We construct the intersection point \(x^*\) and the affine subspace \(E\) in large-scale experiments using a similar method as before, and choose the initial point as \(x_0=P_E(x^*+\alpha G)\), where \(G_{ij}\sim N(0,1)\) are i.i.d. We compare RN-SLRA and Cadzow in nine large-scale experiments; see Fig. 3, where \(n,m\) are the matrix dimensions, \(p=\dim E\), and \(r\) is the target rank. For each tuple \((n,m,p,r)\), both methods use the same initial point, with \(\alpha=0.05\) for the \(100\times100\) cases, \(\alpha=0.02\) for the larger cases, and stopping tolerance \(10^{-8}\).

As shown in Fig. 3, RN-SLRA converges within five iterations in all cases, whereas Cadzow converges much more slowly.

5.3 Low-Rank Approximation of Hankel Matrices↩︎

We now test the performance of RN-SLRA for low-rank approximation of Hankel matrices [11]. We compare RN-SLRA with four representative baseline families for Hankel structured low-rank approximation: Cadzow; VP-SLRA, based on variable projection and reduced nonlinear optimization for affinely structured SLRA [15]; STLN, which approaches the problem from a structured error-correction perspective [10], [46]; and Gradient-system, based on low-rank factor updates with a Hankel penalty [16].

We will construct a single example and then perform both a noise sweep and an outlier sweep to compare RN-SLRA with the baselines on geometrically degenerate Hankel low-rank approximation problems.

5.3.1 Single example↩︎

We consider a \(7\times5\) rank-\(4\) Hankel matrix \(H_c=(\nu_{i+j-1})\), where \(\nu_i=\sum_{\ell=1}^4\beta_\ell z_\ell^i\) for \(i=1,\ldots,11\). Here \(\beta=(1.2,0.9,1.0,1.1)\) and \(z=(e^{-0.2},e^{-0.2+10^{-3}},e^{-0.35},e^{-0.5})\). The closeness of the first two poles makes the instance geometrically degenerate.

The perturbed initial matrix is generated as \[\begin{align} \label{eq:H0} H_0=H_c+\tau \Delta + \eta E_{\text{out}}, \end{align}\tag{40}\] where \(\Delta\) is a random Hankel noise matrix, \(E_{\text{out}}\) is a Hankel matrix supported on the 9-th anti-diagonal, and in this example we choose \(\tau=3\times 10^{-3},\eta=1\times10^{-2}\). \(\Delta\) is a Hankel matrix with generating entries picked independently from \(\mathcal{N}(0,1)\).

We set the parameter of RN-SLRA to \(c=1\times 10^{-7},\rho=1\). We terminate a algorithm if it is stalled or its iteration counts reach the maximal number 300. The numerical results are presented in Fig. 4.

This result shows that RN-SLRA still converges within only a few iterations and rapidly reduces the residual to high accuracy. VP-SLRA also converges to a very small residual, but requires more CPU time than RN-SLRA. Cadzow converges more slowly, while STLN stalls at a moderate residual level, and the gradient-system method fails to make substantial progress within the prescribed maximum number of iterations. These results demonstrate that RN-SLRA remains effective for geometrically degenerate structured low-rank approximation problems, achieving fast convergence and competitive accuracy.

5.3.2 Noise Sweep and Outlier Sweep↩︎

We next test the robustness of RN-SLRA under different noise and outlier levels. In both sweeps, the initial point is generated by 40 , and several independent random Hankel noises \(\Delta\) are used. For the noise sweep, the outlier level \(\eta\) is fixed and \(\tau\in\{10^{-4},3\times10^{-4},10^{-3},3\times10^{-3}\}\). For the outlier sweep, the noise level \(\tau\) is fixed and \(\eta\in\{0,10^{-3},3\times10^{-3},10^{-2}\}\). The mean CPU times over 30 random perturbations are reported in Table 4. Here “fail” means that the algorithm either fails to converge within the prescribed iteration limit or stalls.

4pt

The results show that RN-SLRA remains robust and consistently achieves the smallest mean CPU time across both sweeps. VP-SLRA and Cadzow also converge, but require more time, while STLN and Gradient-system fail on all tested instances. This further illustrates the effectiveness of RN-SLRA on examples with poor geometry.

5.3.3 Large-scale inexact-projection experiment↩︎

We finally compare RN-SLRA with iRN-SLRA when the rank-\(r\) projection is computed only approximately. Here \(n\) and \(m\) denote the matrix dimensions and \(r\) is the target rank. We use larger degenerate Hankel SLRA instances \(H_c^{(n,m,3)}\in\mathbb{R}^{n\times m}\) generated from the repeated-pole model \(\nu_i=\beta_0\lambda^i+\beta_1 i\lambda^i+\beta_2\mu^i\), where the term \(i\lambda^i\) induces poor local geometry. We test \((n,m,r)\in\{(120,160,3),(200,280,3),(350,500,3),(700,1000,3)\}\).

To keep the instances comparable, the poles are moved closer to one as the matrix size increases and the coefficients are mildly rescaled. For these four sizes, the corresponding \((\lambda,\mu,\beta_0,\beta_1,\beta_2)\) are \((e^{-0.06},e^{-0.19},1.0,0.09,0.8)\), \((e^{-0.035},e^{-0.12},1.0,0.06,\) \(0.65)\), \((e^{-0.022},e^{-0.075},1.0,0.045,0.55)\), and \((e^{-0.012},e^{-0.045},1.0,0.025,0.45)\).

The initial matrix is \(H_0=H_c+\tau\Delta+\eta E_{\mathrm{out}}\), where \(\Delta\) is random Hankel noise and \(E_{\mathrm{out}}\) is a Hankel outlier supported on one anti-diagonal. We normalize the perturbations by \(\alpha_c=\|H_c\|/\sqrt{mn}\), and set \(\tau=5\times10^{-2}\alpha_c\), \(\eta=10^{-1}\alpha_c\). For both methods, we tune \(c\in\{10^{-8},3\times10^{-8},10^{-7},3\times10^{-7},10^{-6}\}\) and \(\rho\in\{0,0.5,1.0\}\). For iRN-SLRA, the rank-\(r\) projection is approximated by randomized SVD [38]. We use \(\varepsilon_{\rm r}\in\{0.1,0.25\}\) as the relative tolerance in the randomized SVD routine, and \(s\in\{1,2,3\}\) as the number of subspace iterations. The stopping tolerance is \(10^{-7}\), and the maximum number of outer iterations is \(40\).

Table 5 shows that iRN-SLRA preserves the robustness of RN-SLRA while reducing the running time. The gain becomes more pronounced at larger scales; for \((700,1000,3)\), the time decreases from \(98.14\) to \(58.09\) seconds.

6 Conclusion↩︎

We proposed RN-SLRA, a regularized Newton-type method for local manifold–affine intersection problems. The method is well defined under intrinsic transversality, where we proved local linear convergence, and achieves order \(1+\rho\) under transversality with \(\mu_k=c r_k^\rho\), including quadratic convergence when \(\rho=1\). We also developed an inexact variant, iRN-SLRA, and showed that the same local rates are preserved under a suitable quasioptimal projection condition. Numerical experiments illustrate the robustness of RN-SLRA and the efficiency of iRN-SLRA.

7 Numerical Diagnostics of Intrinsic Transversality↩︎

This appendix reports numerical diagnostics supporting the intrinsic transversality condition in the examples used in the numerical experiments.

7.1 Numerical Principle↩︎

Let \(E\) be the affine constraint set and let \(D_r\) be the fixed-rank manifold of rank-\(r\) matrices. Around a point \(x^*\in E\cap D_r\), we sample nearby points \(x\in E\) and \(y\in D_r\), and define \[u=\frac{x-y}{\|x-y\|_F}.\] We then compute the intrinsic transversality proxy \[\eta(x,y) = \max\Big\{ \operatorname{dist}\big(u,N_{D_r}(y)\big), \operatorname{dist}\big(u,N_E(x)\big) \Big\}.\] Since \(E\) is affine, \(N_E(x)\) is independent of \(x\). A positive lower bound on \(\eta(x,y)\) over nearby sampled pairs indicates that the displacement direction cannot be simultaneously close to both normal geometries, which is the separation expected from intrinsic transversality.

For the small synthetic examples, we also report the singular values of the classical tangent-intersection matrix \(A(x^*)\). A zero or nearly zero smallest singular value indicates failure of classical transversality, whereas positive sampled values of \(\eta\) support the weaker intrinsic transversality condition used in this paper.

7.2 Example 1: Newton-SLRA versus RN-SLRA↩︎

For the first synthetic example, the singular values of the classical tangent-intersection matrix at the solution are \[\sigma(A(x^*)) = (0.9950,\;0.9806,\;0.7645,\;0).\] Thus the classical transversality matrix is rank deficient. Nearby random rank-\(2\) points nevertheless yield positive sampled smallest singular values: \[\begin{array}{c|ccc} \hline \text{sampling scale} & \min \sigma_{\min} & \operatorname{median}\sigma_{\min} & \max \sigma_{\min}\\ \hline 10^{-3} & 5.1222\cdot 10^{-8} & 2.7545\cdot 10^{-6} & 5.0001\cdot 10^{-5}\\ 10^{-2} & 1.2180\cdot 10^{-4} & 3.3858\cdot 10^{-3} & 5.4935\cdot 10^{-2}\\ 5\cdot 10^{-2} & 3.3860\cdot 10^{-5} & 7.0428\cdot 10^{-3} & 1.4893\cdot 10^{-1}\\ \hline \end{array}\]

The intrinsic transversality proxy was then computed for nearby pairs \(x\in E\) and \(y\in D_2\): \[\begin{array}{ccccc} \hline \text{scale}_x & \text{scale}_y & \text{samples} & \min \eta & \operatorname{median}\eta\\ \hline 0.01 & 0.01 & 200 & 0.6611 & 0.9306\\ 0.01 & 0.02 & 200 & 0.8828 & 0.9787\\ 0.01 & 0.05 & 142 & 0.9683 & 0.9937\\ 0.02 & 0.01 & 200 & 0.6706 & 0.8783\\ 0.02 & 0.02 & 200 & 0.7721 & 0.9510\\ 0.02 & 0.05 & 144 & 0.8753 & 0.9896\\ 0.05 & 0.01 & 200 & 0.7323 & 0.9428\\ 0.05 & 0.02 & 200 & 0.7363 & 0.9105\\ 0.05 & 0.05 & 145 & 0.7489 & 0.9585\\ \hline \end{array}\] All sampled minima are bounded away from zero.

7.3 Random Initializations for Example 1↩︎

The random-initialization experiment uses the same affine-rank geometry as Example 1 and changes only the initial point distribution. The intrinsic transversality diagnostics are therefore the same as those reported above.

7.4 Linear-Convergence Example↩︎

For the selected linear-convergence example, the tangent-intersection singular values at \(x^*\) are \[\sigma(A(x^*)) = (0.9968,\;0.9951,\;0.7604,\;9.2694\cdot 10^{-21}),\] so classical transversality again fails numerically. The intrinsic proxy was sampled at scale \(0.02\) with 80 random pairs, giving \[\min \eta = 0.7482.\] This positive value supports intrinsic transversality for this example.

7.5 Cadzow versus RN-SLRA on the Linear-Convergence Example↩︎

The Cadzow comparison uses the same geometry as the linear-convergence example. Hence the corresponding diagnostics are \[\sigma(A(x^*)) = (0.9968,\;0.9951,\;0.7604,\;9.2694\cdot 10^{-21}), \qquad \min \eta = 0.7482.\]

7.6 Large-Scale Synthetic SLRA Examples↩︎

For the nine large-scale synthetic SLRA examples, we report a rank bound for the classical transversality matrix and the minimum sampled intrinsic proxy. The proxy is computed over 24 local samples with sampling radius \(0.01\).

\[\begin{array}{ccccc} \hline (m,n) & p & r & \text{rank bound} & \min \eta\\ \hline (100,100) & 50 & 10 & 2 & 0.9998316\\ (100,100) & 50 & 80 & 2 & 0.9997508\\ (100,100) & 500 & 10 & 2 & 0.9995969\\ (100,100) & 500 & 80 & 2 & 0.9995995\\ (1000,1000) & 50 & 10 & 2 & 0.9999997\\ (1000,1000) & 50 & 100 & 2 & 0.9999983\\ (1000,1000) & 500 & 10 & 2 & 0.9999981\\ (1000,1000) & 500 & 100 & 2 & 0.9999890\\ (1000,100) & 500 & 10 & 2 & 0.9999526\\ \hline \end{array}\]

In all nine cases, the sampled intrinsic proxy is close to one. In contrast, the recorded transversality rank bound is \(2\), much smaller than the number of affine constraints \(p\). This is consistent with the construction: the examples are rank deficient in the classical tangent-intersection sense, while retaining a strong sampled intrinsic separation.

7.7 Hankel SLRA Experiments↩︎

We also computed the intrinsic transversality proxy for the small Hankel experiments. For each case, 200 sampled pairs were used at each radius \(\rho\in\{0.01,0.02,0.05\}\). The table reports the smallest sampled value of \(\eta\) over these radii, together with the range of median values. The column \(\sigma_{\min}(A)\) gives the smallest singular value of the tangent-restriction matrix when it was explicitly formed.

\[\begin{array}{lccccc} \hline \text{experiment} & (m,n) & r & \sigma_{\min}(A) & \min \eta & \operatorname{median}\eta\text{ range}\\ \hline \text{single example} & (7,5) & 4 & 4.5647\cdot 10^{-1} & 0.8628 & [0.9788,0.9800]\\ \text{noise sweep} & (7,5) & 4 & 4.5575\cdot 10^{-1} & 0.8234 & [0.9783,0.9788]\\ \text{noise sweep} & (7,5) & 4 & 4.5647\cdot 10^{-1} & 0.8710 & [0.9788,0.9802]\\ \text{noise sweep} & (8,5) & 2 & 4.8981\cdot 10^{-17} & 0.7718 & [0.8812,0.8838]\\ \text{outlier sweep} & (7,5) & 4 & 4.5575\cdot 10^{-1} & 0.8864 & [0.9803,0.9815]\\ \text{outlier sweep} & (7,5) & 4 & 4.5647\cdot 10^{-1} & 0.8391 & [0.9778,0.9799]\\ \text{outlier sweep} & (8,5) & 2 & 4.8981\cdot 10^{-17} & 0.7627 & [0.8761,0.8792]\\ \hline \end{array}\]

For all small Hankel cases, the sampled intrinsic proxy remains bounded away from zero. Some cases have a nearly singular tangent-restriction matrix, while their intrinsic proxy remains positive. This is consistent with the pattern observed in the synthetic SLRA examples.

7.8 Large-Scale Hankel Intrinsic Diagnostics↩︎

For the large-scale Hankel benchmark used in Table 5, the tangent-restriction matrix was not explicitly formed. The intrinsic proxy was computed directly from sampled normal-cone distances, using 24 sampled pairs at radius \(0.01\) for each case.

\[\begin{array}{ccccc} \hline (m,n) & r & \text{samples} & \min\eta & \operatorname{median}\eta\\ \hline (160,120) & 3 & 24 & 0.9752 & 0.9786\\ (280,200) & 3 & 24 & 0.9852 & 0.9879\\ (500,350) & 3 & 24 & 0.9920 & 0.9927\\ (1000,700) & 3 & 24 & 0.9961 & 0.9964\\ \hline \end{array}\] All four large-scale Hankel cases have sampled \(\eta\)-values close to one, providing further numerical evidence of intrinsic separation.

8 Retraction induced by RN-SLRA in the case \(\rho=0\)↩︎

In this appendix we justify the retraction interpretation in Remark [rmk:retraction-clean-intersection] for the special case \(\rho=0\). In this case the residual-dependent rule \(\mu_k=c r_k^\rho\) is understood as the fixed choice \(\mu_k\equiv c>0\).

Let \(\Phi_c:E\to E\) denote the one-step RN-SLRA map with \(\mu_k\equiv c\). For \(z\in E\) close to \(X\), set \(y=P_M(z)\), \(Q_y:=P_{N_yM}\), and define \[\Phi_c(z) := \arg\min_{u\in E} \left\{ \frac{1}{2}\|u-y\|^2 + \frac{1}{2c}\|Q_y(u-y)\|^2 \right\}.\] We first show that, under clean intersection and suitable smoothness, the limiting map of \(\Phi_c\) defines a genuine retraction on \(X\). We then show that one additional order of smoothness yields a second-order retraction.

References↩︎

---

1. School of Mathematical Sciences, Key Laboratory of the Ministry of Education for Mathematical Foundations and Applications of Digital Technology, University of Science and Technology of China, Hefei, Anhui, China (, ).↩︎

2. Our Python code is available at https://github.com/reniusll/Regularized-Newton-SLRA↩︎