Newton-Schulz Retraction-Based Inference Enables Hidden Quantum Markov Models to Outperform Classical HMMs

# Newton–Schulz Retraction-Based Inference Enables Hidden Quantum Markov Models to Outperform Classical HMMs
Source: [https://arxiv.org/html/2608.06554](https://arxiv.org/html/2608.06554)
Ning Ning

Department of Statistics, Texas A&M University College Station, Texas, USA patning@tamu\.edu

###### Abstract

Hidden Markov models \(HMMs\) are widely used probabilistic models for discrete sequential data but can be limited when hidden dynamics are complex\. Hidden quantum Markov models \(HQMMs\) generalize HMMs by replacing probability vectors with density matrices and stochastic transitions with quantum operations, enabling richer latent representations\. However, existing HQMM learning methods have not consistently outperformed Expectation–Maximization \(EM\)\-trained HMMs on data not generated by quantum processes, limiting their practical applicability\. We introduce NS\-RIS, Newton–Schulz Retraction\-based Inference on the Stiefel manifold, a scalable algorithm for learning trace\-preserving HQMMs\. NS\-RIS uses Newton–Schulz orthogonalization to compute a polar\-factor search direction while preserving Stiefel\-manifold feasibility, avoiding costly matrix decompositions\. We establish a finite\-time stationarity guarantee under standard assumptions on smoothness, stochastic gradients, and finite Newton–Schulz accuracy\. Importantly, NS\-RIS is the first HQMM inference algorithm with mathematical performance guarantee\. Empirically, NS\-RIS provides the first benchmark evidence that an HQMM can significantly outperform an EM\-trained HMM on data not generated by a quantum model\. On synthetic HMM\-generated benchmarks, NS\-RIS outperforms both EM and the state\-of\-the\-art HQMM method COSM, improving the evaluation metric by an average of38\.5%38\.5\\%and by up to50\.6%50\.6\\%\. On a synthetic HQMM benchmark, it improves the test metric over COSM by18\.9%18\.9\\%while reducing runtime by12\.0%12\.0\\%\. On the real\-world Splice classification benchmark, NS\-RIS also surpasses both EM and COSM in higher\-dimensional latent regimes, reducing mean classification error by17\.9%17\.9\\%for latent dimension66and14\.9%14\.9\\%for latent dimension88relative to COSM\. These results move HQMMs beyond a theoretical generalization of HMMs and establish them as practical and expressive models for scientific sequence data\.

###### Contents

1. [1Introduction](https://arxiv.org/html/2608.06554#S1)
2. [2From HMMs to HQMMs](https://arxiv.org/html/2608.06554#S2)1. [2\.1Hidden Markov Models](https://arxiv.org/html/2608.06554#S2.SS1) 2. [2\.2Hidden Quantum Markov Models](https://arxiv.org/html/2608.06554#S2.SS2)
3. [3NS\-RIS for HQMM Learning](https://arxiv.org/html/2608.06554#S3)1. [3\.1The NS\-RIS Algorithm](https://arxiv.org/html/2608.06554#S3.SS1)1. [3\.1\.1Riemannian Incremental Update](https://arxiv.org/html/2608.06554#S3.SS1.SSS1) 2. [3\.1\.2Newton–Schulz Retraction](https://arxiv.org/html/2608.06554#S3.SS1.SSS2) 2. [3\.2Mathematical Guarantees](https://arxiv.org/html/2608.06554#S3.SS2)
4. [4Numerical Analysis](https://arxiv.org/html/2608.06554#S4)1. [4\.1Setup and Evaluation Metric](https://arxiv.org/html/2608.06554#S4.SS1) 2. [4\.2Synthetic HMM Test Performance](https://arxiv.org/html/2608.06554#S4.SS2) 3. [4\.3Synthetic HQMM Test Performance](https://arxiv.org/html/2608.06554#S4.SS3)
5. [5Empirical Analysis](https://arxiv.org/html/2608.06554#S5)1. [5\.1The Splice Dataset](https://arxiv.org/html/2608.06554#S5.SS1) 2. [5\.2Splice Classification Results](https://arxiv.org/html/2608.06554#S5.SS2) 3. [5\.3Scientific Validation](https://arxiv.org/html/2608.06554#S5.SS3)
6. [6Conclusion](https://arxiv.org/html/2608.06554#S6)
7. [AProofs](https://arxiv.org/html/2608.06554#A1)1. [A\.1Proof of Proposition1](https://arxiv.org/html/2608.06554#A1.SS1) 2. [A\.2Proof of Theorem1](https://arxiv.org/html/2608.06554#A1.SS2)
8. [References](https://arxiv.org/html/2608.06554#bib)

## 1Introduction

Hidden Markov models \(HMMs\) are a standard language for sequential data\. They are used wherever an observed time series is driven by an unobserved state process, including computational biology, speech and language processing, signal analysis, finance, and many other scientific areas\. In biology alone, HMMs and profile HMMs have become foundational tools for modeling protein families, sequence motifs, and splice\-site structure\(Krogh et al\.,,[1994](https://arxiv.org/html/2608.06554#bib.bib17); Eddy,,[1998](https://arxiv.org/html/2608.06554#bib.bib9); Burge and Karlin,,[1997](https://arxiv.org/html/2608.06554#bib.bib5)\)\. Their success comes from a simple and interpretable architecture: a hidden Markov chain evolves over time and emits observations conditionally on its current state\. This simplicity also creates a limitation\. When the hidden mechanism contains higher\-order, contextual, or nonclassical dependencies, a probability vector over a finite set of latent states may require a large state space or may fail to represent the relevant dependence efficiently\.

Hidden quantum Markov models \(HQMMs\) provide a principled extension of this classical framework\. More broadly, researchers in physics and machine learning have developed quantum graphical models by incorporating the quantum mechanical view of probability into graphical\-model inference\(Warmuth and Kuzmin,,[2006](https://arxiv.org/html/2608.06554#bib.bib33); Leifer and Poulin,,[2008](https://arxiv.org/html/2608.06554#bib.bib18); Yeang,,[2010](https://arxiv.org/html/2608.06554#bib.bib34); Leifer and Spekkens,,[2013](https://arxiv.org/html/2608.06554#bib.bib19)\)\. Instead of representing the latent belief state by a probability vector, an HQMM represents it by a density matrix; instead of using nonnegative transition\-emission matrices, it uses symbol\-conditioned quantum operations represented by Kraus operators\(Monras et al\.,,[2010](https://arxiv.org/html/2608.06554#bib.bib20); Clark et al\.,,[2015](https://arxiv.org/html/2608.06554#bib.bib6); Deb et al\.,,[2026](https://arxiv.org/html/2608.06554#bib.bib7)\)\. This formulation preserves the sequential likelihood structure of HMMs while allowing the latent state to encode coherence through off\-diagonal matrix entries\. Prior work has shown that HQMMs can be more expressive than classical HMMs and can represent some sequential processes more compactly\([Srinivasan et al\., 2018a,](https://arxiv.org/html/2608.06554#bib.bib29); Ning,,[2025](https://arxiv.org/html/2608.06554#bib.bib24)\)\. Recent applications and physical formulations have further connected HQMMs to sequential analysis, measurement\-induced quantum inference, and many\-body topological structure\(Souissi and Andolsi,,[2026](https://arxiv.org/html/2608.06554#bib.bib28); Kim et al\.,,[2026](https://arxiv.org/html/2608.06554#bib.bib15)\)\. General overviews of quantum machine learning also identify learning expressive quantum sequential models from data as an important open problem\(Schuld et al\.,,[2015](https://arxiv.org/html/2608.06554#bib.bib27); Biamonte et al\.,,[2016](https://arxiv.org/html/2608.06554#bib.bib4)\)\. Related quantum graphical model perspectives connect HQMMs to inference in Hilbert space and to operator\-based models of stochastic processes\(Jaeger,,[2000](https://arxiv.org/html/2608.06554#bib.bib13); Zhu et al\.,,[2025](https://arxiv.org/html/2608.06554#bib.bib35)\)\.

The central obstacle is learning\. HQMM parameters must satisfy the trace\-preserving Kraus constraint∑y,qKy,q†​Ky,q=I\\sum\_\{y,q\}K\_\{y,q\}^\{\\dagger\}K\_\{y,q\}=I, which makes the feasible parameter space a complex Stiefel manifold after stacking the Kraus operators\. Existing methods either pay a high computational price to maintain feasibility or use updates that do not consistently translate HQMM expressiveness into better performance on ordinary, non\-quantum\-generated data\. As a result, HQMMs have remained more compelling as a theoretical generalization of HMMs than as a practical replacement for HMMs in the broad scientific settings where HMMs are routinely used\.

This paper introduces NS\-RIS, short for Newton–Schulz Retraction\-Based Inference on the Stiefel manifold, for scalable HQMM learning\. The method uses Newton–Schulz orthogonalization twice: first to approximate the polar factor of the Riemannian momentum direction, and second to retract the updated Kraus matrix back to the Stiefel manifold\. In this way, NS\-RIS preserves the physical trace\-preserving constraint while avoiding expensive matrix decompositions inside the repeated learning loop\. We compare NS\-RIS with the main existing HQMM learning baselines: Givens Search \(GS\)\([Srinivasan et al\., 2018b,](https://arxiv.org/html/2608.06554#bib.bib30)\), constrained optimization on the Stiefel manifold \(COSM\)\(Adhikary et al\.,,[2020](https://arxiv.org/html/2608.06554#bib.bib1)\), and the classical Expectation–Maximization \(EM\) procedure for HMMs which serves as the standard non\-quantum benchmark\.

The key findings are as follows\.

- •NS\-RIS gives a decomposition\-free learning method on the Stiefel manifold for trace\-preserving HQMMs\. Algorithm[1](https://arxiv.org/html/2608.06554#alg1)gives the full incremental learning procedure and Algorithm[2](https://arxiv.org/html/2608.06554#alg2)gives the Newton–Schulz orthogonalization subroutine; together they make the Kraus constraint part of the geometry of the algorithm rather than an after\-the\-fact correction\.
- •Proposition[1](https://arxiv.org/html/2608.06554#Thmproposition1)mathematically explains the two purposes of the Newton–Schulz design\. The first Newton–Schulz step approximates the polar factor that solves the operator\-norm steepest\-descent subproblem under constraint in the \(matrix\) operator norm, with optimal value obtained with the \(matrix\) nuclear norm; the second Newton–Schulz step approximates the nearest Stiefel retraction under the \(matrix\) Frobenius norm\.
- •The main theoretical contribution is Theorem[1](https://arxiv.org/html/2608.06554#Thmtheorem1), which establishes a finite\-time stationarity guarantee under standard smoothness, stochastic\-gradient, and finite Newton–Schulz accuracy assumptions\. This appears to be the first mathematical convergence guarantee for an HQMM inference algorithm\. The bound separates the descent term, smoothness term, finite Newton–Schulz feasibility residuals, and stochastic momentum\-tracking term, showing how each source affects the final stationarity level\.
- •On synthetic HMM\-generated data, generated from the standard66\-hidden\-state,66\-output HMM benchmark, NS\-RIS exceeds the EM\-trained HMM baseline and improves over COSM by38\.5%38\.5\\%on average across tested latent dimensions and Kraus ranks\. In the best\-performing configuration, it reaches a50\.6%50\.6\\%relative improvement over COSM while maintaining competitive running time; see Subsection[4\.2](https://arxiv.org/html/2608.06554#S4.SS2)and Figures[1](https://arxiv.org/html/2608.06554#S4.F1)–[2](https://arxiv.org/html/2608.06554#S4.F2)\.
- •On the synthetic HQMM benchmark generated by a22\-hidden\-state,66\-output HQMM with Kraus rankw=1w=1, NS\-RIS achieves the best training, validation, and test metrics among GS, COSM, and EM\. It improves the test metric over COSM by18\.9%18\.9\\%while reducing runtime by12\.0%12\.0\\%, showing that the gain is not obtained by a larger computational cost; see Subsection[4\.3](https://arxiv.org/html/2608.06554#S4.SS3)and Figure[3](https://arxiv.org/html/2608.06554#S4.F3)\.
- •On the real Splice classification benchmark, NS\-RIS outperforms EM and COSM in higher\-dimensional latent regimes, reducing mean error by17\.9%17\.9\\%for latent dimension66and by14\.9%14\.9\\%for latent dimension88relative to COSM after averaging over Kraus ranks\. Class\-wise results further show improvements in EI, IE, and negative\-class errors for\(n,w\)=\(6,4\)\(n,w\)=\(6,4\)and\(8,4\)\(8,4\), with the largest reductions in the negative\-class error\. The scientific validation in Subsection[5\.3](https://arxiv.org/html/2608.06554#S5.SS3)explains why latent dimensions larger than the four\-letter nucleotide alphabet are biologically meaningful: they can encode splice motifs, exon–intron regimes, positional context, and long\-range dependencies; see Figures[4](https://arxiv.org/html/2608.06554#S5.F4)–[6](https://arxiv.org/html/2608.06554#S5.F6)\.

Together, these findings provide the first benchmark evidence in this setting that an HQMM can substantially outperform an EM\-trained HMM on data not generated by a quantum model; they also establish the first HQMM inference algorithm with a mathematical convergence guarantee\. This moves HQMMs beyond the statement that they generalize HMMs mathematically: it shows that, with an effective Stiefel\-manifold learning procedure, HQMMs can reliably deliver practical gains on ordinary scientific sequence data\.

The remainder of the paper is organized as follows\. Section[2](https://arxiv.org/html/2608.06554#S2)reviews HMMs and introduces HQMMs in the operator form used throughout the paper\. Section[3](https://arxiv.org/html/2608.06554#S3)formulates HQMM learning on the Stiefel manifold, presents NS\-RIS, and states the main mathematical guarantees\. Section[4](https://arxiv.org/html/2608.06554#S4)evaluates NS\-RIS on synthetic HMM and HQMM benchmarks\. Section[5](https://arxiv.org/html/2608.06554#S5)studies real splice sequence classification and discusses its scientific implications\. Section[6](https://arxiv.org/html/2608.06554#S6)concludes, and Appendix[A](https://arxiv.org/html/2608.06554#A1)contains the proofs\.

## 2From HMMs to HQMMs

This section recalls the classical HMM formulation in Subsection[2\.1](https://arxiv.org/html/2608.06554#S2.SS1), and then introduces the HQMM extension in Subsection[2\.2](https://arxiv.org/html/2608.06554#S2.SS2)whose quantum states and Kraus operators generalize the hidden\-state dynamics\.

### 2\.1Hidden Markov Models

An HMM models a discrete\-time sequence by separating what is observed from what is only indirectly inferred\. At timett, the model has a hidden stateXt∈\{1,…,n\}X\_\{t\}\\in\\\{1,\\ldots,n\\\}and emits an observationYt∈\{1,…,m\}Y\_\{t\}\\in\\\{1,\\ldots,m\\\}\. The hidden process is Markovian, so the distribution ofXtX\_\{t\}depends on the past only throughXt−1X\_\{t\-1\}, while the observationYtY\_\{t\}is drawn conditionally on the current hidden stateXtX\_\{t\}\. This architecture is useful when the data exhibit temporal dependence but the mechanism driving that dependence is not directly observed\.

Letπ∈ℝn\\pi\\in\\mathbb\{R\}^\{n\}be the initial distribution, letA∈ℝ≥0n×nA\\in\\mathbb\{R\}\_\{\\geq 0\}^\{n\\times n\}be the transition matrix withAi​j=P​\(Xt=i∣Xt−1=j\)A\_\{ij\}=P\(X\_\{t\}=i\\mid X\_\{t\-1\}=j\), and letC∈ℝ≥0m×nC\\in\\mathbb\{R\}\_\{\\geq 0\}^\{m\\times n\}be the emission matrix withCy​i=P​\(Yt=y∣Xt=i\)C\_\{yi\}=P\(Y\_\{t\}=y\\mid X\_\{t\}=i\)\. We use the column\-stochastic convention𝟏T​π=1\\mathbf\{1\}^\{T\}\\pi=1,𝟏T​A=𝟏T\\mathbf\{1\}^\{T\}A=\\mathbf\{1\}^\{T\}, and𝟏T​C=𝟏T\\mathbf\{1\}^\{T\}C=\\mathbf\{1\}^\{T\}\. Ifxt−1x\_\{t\-1\}denotes the filtering distribution over hidden states after the firstt−1t\-1observations, prediction givesA​xt−1Ax\_\{t\-1\}before the next symbol is observed\. OnceYt=ytY\_\{t\}=y\_\{t\}is observed, the filtering distribution is updated by

xt=diag⁡\(C​\(yt,:\)\)​A​xt−1𝟏T​diag⁡\(C​\(yt,:\)\)​A​xt−1\.x\_\{t\}=\\frac\{\\operatorname\{diag\}\(C\(y\_\{t\},:\)\)Ax\_\{t\-1\}\}\{\\mathbf\{1\}^\{T\}\\operatorname\{diag\}\(C\(y\_\{t\},:\)\)Ax\_\{t\-1\}\}\.The denominator is the one\-step predictive probability ofyty\_\{t\}under the current model\.

It is often convenient to combine transition and emission into a symbol\-indexed operator

Ty=diag⁡\(C​\(y,:\)\)​A\.T\_\{y\}=\\operatorname\{diag\}\(C\(y,:\)\)A\.\(1\)For an observation sequencey¯=y1,…,yT\\bar\{y\}=y\_\{1\},\\ldots,y\_\{T\}, the likelihood is then

P​\(y¯\)=𝟏T​TyT​TyT−1​⋯​Ty1​π\.P\(\\bar\{y\}\)=\\mathbf\{1\}^\{T\}T\_\{y\_\{T\}\}T\_\{y\_\{T\-1\}\}\\cdots T\_\{y\_\{1\}\}\\pi\.Thus an HMM can be viewed as a family of nonnegative linear maps\{Ty\}y=1m\\\{T\_\{y\}\\\}\_\{y=1\}^\{m\}acting on probability vectors, with normalization after each observation when filtering is required\. This operator view is the form that most directly connects classical HMMs to their quantum generalization: HQMMs keep the same sequence\-likelihood logic, but replace probability vectors by density matrices and replace the nonnegative mapsTyT\_\{y\}by quantum operations\.

![[Uncaptioned image]](https://arxiv.org/html/2608.06554v1/pic/HMM_HQMM.png)
### 2\.2Hidden Quantum Markov Models

The passage from HMMs to HQMMs begins with the state representation\. In a classical model, the belief state is a probability vector\. In a quantum model, the state is represented by a density matrixρ\\rho, a positive semidefinite matrix withtr⁡\(ρ\)=1\\operatorname\{tr\}\(\\rho\)=1\. The diagonal entries ofρ\\rhocan be interpreted as ordinary probabilities in a chosen basis, while the off\-diagonal entries encode quantum coherence\. For example, the two\-dimensional state

\|ψ⟩=12​\|0⟩−i2​\|1⟩\|\\psi\\rangle=\\frac\{1\}\{\\sqrt\{2\}\}\|0\\rangle\-\\frac\{i\}\{\\sqrt\{2\}\}\|1\\ranglehas density matrix

ρ=\|ψ⟩​⟨ψ\|=\[1/2i/2−i/21/2\]\.\\rho=\|\\psi\\rangle\\langle\\psi\|=\\begin\{bmatrix\}1/2&i/2\\\\ \-i/2&1/2\\end\{bmatrix\}\.The diagonal entries assign equal probability to the two basis states, and the off\-diagonal entries retain phase information that is absent from a classical probability vector\.

The second replacement is at the level of dynamics\. A quantum operation𝒦\\mathcal\{K\}maps density matrices to density matrices and can be written in Kraus form as

𝒦​\(ρ\)=∑iKi​ρ​Ki†,\\mathcal\{K\}\(\\rho\)=\\sum\_\{i\}K\_\{i\}\\rho K\_\{i\}^\{\\dagger\},whereKi†K\_\{i\}^\{\\dagger\}is the Hermitian conjugate ofKiK\_\{i\}\. The operation is trace\-preserving when

∑iKi†​Ki=𝕀n,\\sum\_\{i\}K\_\{i\}^\{\\dagger\}K\_\{i\}=\\mathbb\{I\}\_\{n\},and trace\-nonincreasing when the left\-hand side is bounded above by𝕀n\\mathbb\{I\}\_\{n\}\. In an HQMM, each observable symbol is associated with a trace\-nonincreasing operation, and the sum over all symbols is trace\-preserving\. This mirrors the HMM requirement that the probabilities of all possible emissions sum to one\.

###### Definition 1\(Monras et al\., \([2010](https://arxiv.org/html/2608.06554#bib.bib20)\)\)\.

A hidden quantum Markov model is a quantum system with stateρ\\rhoand a set of quantum operations\{𝒦s\}\\\{\\mathcal\{K\}\_\{s\}\\\}indexed by output symbols, such that∑s𝒦s\\sum\_\{s\}\\mathcal\{K\}\_\{s\}is trace\-preserving\. At each time step, symbolssis generated with probabilityP​\(s\)=tr⁡\[𝒦s​\(ρ\)\]P\(s\)=\\operatorname\{tr\}\[\\mathcal\{K\}\_\{s\}\(\\rho\)\], and the state is updated toρs=𝒦s​\(ρ\)/P​\(s\)\\rho\_\{s\}=\\mathcal\{K\}\_\{s\}\(\\rho\)/P\(s\)\.

If each symbol\-conditioned operation is represented bywwKraus operators,

𝒦s​\(ρ\)=∑q=1wKs,q​ρ​Ks,q†,\\mathcal\{K\}\_\{s\}\(\\rho\)=\\sum\_\{q=1\}^\{w\}K\_\{s,q\}\\rho K\_\{s,q\}^\{\\dagger\},then the trace\-preserving condition for the full model is

∑s=1m∑q=1wKs,q†​Ks,q=𝕀n\.\\sum\_\{s=1\}^\{m\}\\sum\_\{q=1\}^\{w\}K\_\{s,q\}^\{\\dagger\}K\_\{s,q\}=\\mathbb\{I\}\_\{n\}\.\(2\)Given an observed sequencey1,…,yTy\_\{1\},\\ldots,y\_\{T\}, the HQMM filtering recursion is

ρt=∑q=1wKyt,q​ρt−1​Kyt,q†tr⁡\(∑q=1wKyt,q​ρt−1​Kyt,q†\)\.\\rho\_\{t\}=\\frac\{\\sum\_\{q=1\}^\{w\}K\_\{y\_\{t\},q\}\\rho\_\{t\-1\}K\_\{y\_\{t\},q\}^\{\\dagger\}\}\{\\operatorname\{tr\}\\left\(\\sum\_\{q=1\}^\{w\}K\_\{y\_\{t\},q\}\\rho\_\{t\-1\}K\_\{y\_\{t\},q\}^\{\\dagger\}\\right\)\}\.\(3\)The denominator is the predictive probability of the emitted symbol\. Thus HQMMs retain the forward\-propagation structure of HMMs while allowing the latent state to carry quantum coherence and allowing the transition\-emission mechanism to be modeled by completely positive maps\. This additional structure is one reason HQMMs can represent some sequential dependencies more compactly than classical HMMs\([Srinivasan et al\., 2018b,](https://arxiv.org/html/2608.06554#bib.bib30); Adhikary et al\.,,[2020](https://arxiv.org/html/2608.06554#bib.bib1)\)\.

## 3NS\-RIS for HQMM Learning

This section first formulates HQMM learning on the Stiefel manifold and describes the NS\-RIS algorithm in Subsection[3\.1](https://arxiv.org/html/2608.06554#S3.SS1)\. Subsection[3\.2](https://arxiv.org/html/2608.06554#S3.SS2)then summarizes the main mathematical assumptions and convergence guarantees; the detailed proofs are deferred to Appendix[A](https://arxiv.org/html/2608.06554#A1)\.

### 3\.1The NS\-RIS Algorithm

Learning an HQMM requires estimating a collection of symbol\-conditioned Kraus operators that explain the observed sequences while satisfying the trace\-preserving constraint that makes the model physically admissible\. LetDtr=\{Bi\}i=1ND\_\{\\mathrm\{tr\}\}=\\\{B\_\{i\}\\\}\_\{i=1\}^\{N\}denote the training set, whereBi=\(yi​1,…,yi​Ti\)B\_\{i\}=\(y\_\{i1\},\\ldots,y\_\{iT\_\{i\}\}\)is a discrete observation sequence\. For a single sequenceB=\(y1,…,yT\)B=\(y\_\{1\},\\ldots,y\_\{T\}\)and initial density matrixρ0\\rho\_\{0\}, the log\-likelihood is obtained by propagating the unnormalized quantum state through the Kraus operators associated with the observed symbols:

ℓ​\(K;B\)=log⁡tr⁡\(∑qT=1wKyT,qT​⋯​\(∑q1=1wKy1,q1​ρ0​Ky1,q1†\)​⋯​KyT,qT†\)\.\\ell\(K;B\)=\\log\\operatorname\{tr\}\\left\(\\sum\_\{q\_\{T\}=1\}^\{w\}K\_\{y\_\{T\},q\_\{T\}\}\\cdots\\left\(\\sum\_\{q\_\{1\}=1\}^\{w\}K\_\{y\_\{1\},q\_\{1\}\}\\rho\_\{0\}K\_\{y\_\{1\},q\_\{1\}\}^\{\\dagger\}\\right\)\\cdots K\_\{y\_\{T\},q\_\{T\}\}^\{\\dagger\}\\right\)\.\(4\)Here,KKdenotes the full Kraus tensor\{Ky,q:y=1,…,s,q=1,…,w\}\\\{K\_\{y,q\}:y=1,\\ldots,s,\\ q=1,\\ldots,w\\\}, wheressis the output alphabet size,wwis the number of Kraus operators per output, and eachKy,q∈ℂn×nK\_\{y,q\}\\in\\mathbb\{C\}^\{n\\times n\}acts on ann×nn\\times n\-dimensional latent quantum state\. The learning objective is the empirical negative log\-likelihood

F​\(K\)=1N​∑i=1NL​\(K;Bi\),L​\(K;Bi\)=−ℓ​\(K;Bi\)\.F\(K\)=\\frac\{1\}\{N\}\\sum\_\{i=1\}^\{N\}L\(K;B\_\{i\}\),\\qquad L\(K;B\_\{i\}\)=\-\\ell\(K;B\_\{i\}\)\.
Direct optimization ofFFis challenging for two reasons\. First, the likelihood in \([4](https://arxiv.org/html/2608.06554#S3.E4)\) involves repeated products of noncommuting complex matrices\. Second, the Kraus operators must jointly satisfy

∑y=1s∑q=1wKy,q†​Ky,q=𝕀n,\\sum\_\{y=1\}^\{s\}\\sum\_\{q=1\}^\{w\}K\_\{y,q\}^\{\\dagger\}K\_\{y,q\}=\\mathbb\{I\}\_\{n\},which is a nonlinear matrix constraint\. NS\-RIS handles this constraint geometrically by stacking all Kraus operators into a single matrixΓ=\(K1,1T,…,Ks,wT\)T∈ℂs​w​n×n\\Gamma=\(K\_\{1,1\}^\{T\},\\ldots,K\_\{s,w\}^\{T\}\)^\{T\}\\in\\mathbb\{C\}^\{swn\\times n\}\. Withp=s​w​np=swn, the trace\-preserving condition is equivalent toΓ†​Γ=𝕀n\\Gamma^\{\\dagger\}\\Gamma=\\mathbb\{I\}\_\{n\}, so the feasible parameter space is the complex Stiefel manifold

St​\(p,n\)=\{Γ∈ℂp×n:Γ†​Γ=𝕀n\}\.\\mathrm\{St\}\(p,n\)=\\\{\\Gamma\\in\\mathbb\{C\}^\{p\\times n\}:\\Gamma^\{\\dagger\}\\Gamma=\\mathbb\{I\}\_\{n\}\\\}\.Thus HQMM learning can be formulated as constrained optimization overSt​\(p,n\)\\mathrm\{St\}\(p,n\)\.

#### 3\.1\.1Riemannian Incremental Update

LetΓk\\Gamma\_\{k\}denote the stacked Kraus matrix at thekk\-th inner iteration\. NS\-RIS selects a mini\-batch index setℬk\\mathcal\{B\}\_\{k\}according to the current random permutation of the training set; the incremental single\-sequence case corresponds to cardinality\|ℬk\|=1\|\\mathcal\{B\}\_\{k\}\|=1\. For mathematical convenience, we write the sequence loss asL​\(Γk;Bik\)L\(\\Gamma\_\{k\};B\_\{i\_\{k\}\}\)forik∈ℬki\_\{k\}\\in\\mathcal\{B\}\_\{k\}\. This denotes the same negative log\-likelihoodL​\(Kk;Bik\)=−ℓ​\(Kk;Bik\)L\(K\_\{k\};B\_\{i\_\{k\}\}\)=\-\\ell\(K\_\{k\};B\_\{i\_\{k\}\}\)defined above, after identifying the Kraus tensorKkK\_\{k\}with its stacked Stiefel representationΓk\\Gamma\_\{k\}\. Thus, wheneverL​\(Γk;Bik\)L\(\\Gamma\_\{k\};B\_\{i\_\{k\}\}\)appears, the Kraus operators used in the likelihood \([4](https://arxiv.org/html/2608.06554#S3.E4)\) are the block components ofΓk\\Gamma\_\{k\}\. The Euclidean gradient used at this iteration is

Gk=1\|ℬk\|​∑ik∈ℬk∇ΓL​\(Γk;Bik\)∈ℂp×n\.G\_\{k\}=\\frac\{1\}\{\|\\mathcal\{B\}\_\{k\}\|\}\\sum\_\{i\_\{k\}\\in\\mathcal\{B\}\_\{k\}\}\\nabla\_\{\\Gamma\}L\(\\Gamma\_\{k\};B\_\{i\_\{k\}\}\)\\in\\mathbb\{C\}^\{p\\times n\}\.SinceΓk\\Gamma\_\{k\}is constrained to the complex Stiefel manifold, the update direction must lie in the tangent space

TΓk​St​\(p,n\)=\{ξ∈ℂp×n:Γk†​ξ\+ξ†​Γk=0\};T\_\{\\Gamma\_\{k\}\}\\mathrm\{St\}\(p,n\)=\\\{\\xi\\in\\mathbb\{C\}^\{p\\times n\}:\\Gamma\_\{k\}^\{\\dagger\}\\xi\+\\xi^\{\\dagger\}\\Gamma\_\{k\}=0\\\};This tangent\-space equation is the first\-order linearization of the Stiefel constraintΓk†​Γk=𝕀n\\Gamma\_\{k\}^\{\\dagger\}\\Gamma\_\{k\}=\\mathbb\{I\}\_\{n\}\. NS\-RIS uses the orthogonal projection

ΠΓk​\(Gk\)=Gk−12​Γk​\(Gk†​Γk\+Γk†​Gk\),\\Pi\_\{\\Gamma\_\{k\}\}\(G\_\{k\}\)=G\_\{k\}\-\\frac\{1\}\{2\}\\Gamma\_\{k\}\\left\(G\_\{k\}^\{\\dagger\}\\Gamma\_\{k\}\+\\Gamma\_\{k\}^\{\\dagger\}G\_\{k\}\\right\),\(5\)which belongs toTΓk​St​\(p,n\)T\_\{\\Gamma\_\{k\}\}\\mathrm\{St\}\(p,n\)wheneverΓk†​Γk=𝕀n\\Gamma\_\{k\}^\{\\dagger\}\\Gamma\_\{k\}=\\mathbb\{I\}\_\{n\}\. Rather than forming a full gradient over all training sequences, NS\-RIS processes sequences incrementally and maintains a momentum\-smoothed tangent direction

Mk\+1=β​Mk\+\(1−β\)​ΠΓk​\(Gk\),0≤β<1\.M\_\{k\+1\}=\\beta M\_\{k\}\+\(1\-\\beta\)\\Pi\_\{\\Gamma\_\{k\}\}\(G\_\{k\}\),\\qquad 0\\leq\\beta<1\.
![[Uncaptioned image]](https://arxiv.org/html/2608.06554v1/pic/algorithm.png)
#### 3\.1\.2Newton–Schulz Retraction

The retraction step is implemented using a Newton–Schulz approximate orthogonalization map, denoted byNSOrth⁡\(⋅,TNS\)\\operatorname\{NSOrth\}\(\\cdot,T\_\{\\mathrm\{NS\}\}\), whereTNST\_\{\\mathrm\{NS\}\}is the prescribed number of Newton–Schulz iterations\. Given a matrixZ∈ℂp×nZ\\in\\mathbb\{C\}^\{p\\times n\}, the procedure first rescalesZ0=Z‖Z‖F\+εZ\_\{0\}=\\frac\{Z\}\{\\\|Z\\\|\_\{F\}\+\\varepsilon\}, where∥⋅∥F\\\|\\cdot\\\|\_\{F\}denotes the Frobenius norm, and then applies, fort=0,…,TNS−1t=0,\\ldots,T\_\{\\mathrm\{NS\}\}\-1,

Ht=Zt†​Zt,Zt\+1=a​Zt\+Zt​\(b​Ht\+c​Ht2\),H\_\{t\}=Z\_\{t\}^\{\\dagger\}Z\_\{t\},\\qquad Z\_\{t\+1\}=aZ\_\{t\}\+Z\_\{t\}\(bH\_\{t\}\+cH\_\{t\}^\{2\}\),with fixed coefficientsa=3\.4445a=3\.4445,b=−4\.7750b=\-4\.7750, andc=2\.0315c=2\.0315\. The outputNSOrth⁡\(Z,TNS\)=ZTNS\\operatorname\{NSOrth\}\(Z,T\_\{\\mathrm\{NS\}\}\)=Z\_\{T\_\{\\mathrm\{NS\}\}\}serves as a decomposition\-free approximate retraction toward the Stiefel manifold\.

At iterationkk, NS\-RIS uses this map in two distinct ways; rigorous justification is provided in Proposition[1](https://arxiv.org/html/2608.06554#Thmproposition1)\. First, it is applied to the momentum buffer to obtain a polar\-like, normalized search direction,

M~k\+1=NSOrth⁡\(Mk\+1,TNS\),\\widetilde\{M\}\_\{k\+1\}=\\operatorname\{NSOrth\}\(M\_\{k\+1\},T\_\{\\mathrm\{NS\}\}\),which plays the role of an orthogonalized descent direction\. This step is analogous to replacing the raw momentum by its polar factor, but avoids an explicit singular value decomposition\. Second, after moving along this direction, the intermediate matrix generally no longer satisfies the Stiefel constraint\. Indeed, even ifΓk†​Γk=𝕀n\\Gamma\_\{k\}^\{\\dagger\}\\Gamma\_\{k\}=\\mathbb\{I\}\_\{n\}andM~k\+1\\widetilde\{M\}\_\{k\+1\}has approximately orthonormal columns, the matrixΓk\+12=Γk−ηk​M~k\+1\\Gamma\_\{k\+\\frac\{1\}\{2\}\}=\\Gamma\_\{k\}\-\\eta\_\{k\}\\widetilde\{M\}\_\{k\+1\}need not satisfyΓk\+12†​Γk\+12=𝕀n\\Gamma\_\{k\+\\frac\{1\}\{2\}\}^\{\\dagger\}\\Gamma\_\{k\+\\frac\{1\}\{2\}\}=\\mathbb\{I\}\_\{n\}\. Therefore, NS\-RIS applies the same Newton–Schulz orthogonalization a second time as an efficient approximation to the polar projection of the updated parameter matrix back to the feasible Stiefel manifold:

Γk\+12=Γk−ηk​M~k\+1,Γk\+1=NSOrth⁡\(Γk\+12,TNS\)\.\\Gamma\_\{k\+\\frac\{1\}\{2\}\}=\\Gamma\_\{k\}\-\\eta\_\{k\}\\widetilde\{M\}\_\{k\+1\},\\qquad\\Gamma\_\{k\+1\}=\\operatorname\{NSOrth\}\(\\Gamma\_\{k\+\\frac\{1\}\{2\}\},T\_\{\\mathrm\{NS\}\}\)\.Thus, the first call toNSOrth\\operatorname\{NSOrth\}controls the geometry of the update direction, whereas the second call enforces the constraintΓk\+1†​Γk\+1=𝕀n\\Gamma\_\{k\+1\}^\{\\dagger\}\\Gamma\_\{k\+1\}=\\mathbb\{I\}\_\{n\}required by the stacked Kraus representation of a valid HQMM\.

The updated matrixΓk\+1\\Gamma\_\{k\+1\}is finally reshaped back into Kraus\-operator form\. At the end of each epoch, the stepsize is decayed geometrically,ηe\+1=α​ηe\\eta\_\{e\+1\}=\\alpha\\eta\_\{e\}, for0<α≤10<\\alpha\\leq 1, whereα\\alphais a fixed decay factor controlling how quickly the learning rate decreases across epochs\. When validation data are available, the parameter tensor with the largest validation log\-likelihood is retained:

Kbest=arg⁡maxKe⁡ℓval​\(Ke\)\.K\_\{\\mathrm\{best\}\}=\\arg\\max\_\{K\_\{e\}\}\\ell\_\{\\mathrm\{val\}\}\(K\_\{e\}\)\.The complete NS\-RIS training procedure is summarized in Algorithm[1](https://arxiv.org/html/2608.06554#alg1)\. The Newton–Schulz orthogonalization subroutine used in Algorithm[1](https://arxiv.org/html/2608.06554#alg1)is given explicitly in Algorithm[2](https://arxiv.org/html/2608.06554#alg2)\.

Algorithm 1NS\-RIS for HQMM Learning1:Training sequences

DtrD\_\{\\mathrm\{tr\}\}, initial Kraus tensor

K0K\_\{0\}, initial stepsize

η0\\eta\_\{0\}, momentum parameter

β\\beta, decay factor

α\\alpha, Newton–Schulz iteration count

TNST\_\{\\mathrm\{NS\}\}, number of epochs

EE, optional validation set

DvalD\_\{\\mathrm\{val\}\}
2:Best Kraus tensor

KbestK\_\{\\mathrm\{best\}\}
3:Convert

K0K\_\{0\}to a Stiefel matrix

Γ∈St​\(p,n\)\\Gamma\\in\\mathrm\{St\}\(p,n\), where

p=s​w​np=swn
4:Set

Kbest←K0K\_\{\\mathrm\{best\}\}\\leftarrow K\_\{0\},

ℓbest←−∞\\ell\_\{\\mathrm\{best\}\}\\leftarrow\-\\infty
5:Set stepsize

η←η0\\eta\\leftarrow\\eta\_\{0\}
6:Initialize momentum buffer

M←0p×nM\\leftarrow 0\_\{p\\times n\}
7:Record initial train and validation metrics in

ℋ\\mathcal\{H\}
8:for

e=1,…,Ee=1,\\ldots,Edo

9:Randomly permute the training sequence indices

10:Reset epoch running totals

11:foreach sequence index

iiin the permutationdo

12:Let

B←Dtr​\[i,:\]B\\leftarrow D\_\{\\mathrm\{tr\}\}\[i,:\]
13:Compute the gradient tensor

∇KL​\(K;B\)\\nabla\_\{K\}L\(K;B\)using the HQMM gradient routine

14:Convert

∇KL​\(K;B\)\\nabla\_\{K\}L\(K;B\)to a matrix

G∈ℂp×nG\\in\\mathbb\{C\}^\{p\\times n\}
15:Project

GGto the tangent space at

Γ\\Gamma:

GR=G−12​Γ​\(G†​Γ\+Γ†​G\)G\_\{R\}=G\-\\frac\{1\}\{2\}\\Gamma\\left\(G^\{\\dagger\}\\Gamma\+\\Gamma^\{\\dagger\}G\\right\)
16:Update the momentum buffer:

M←β​M\+\(1−β\)​GRM\\leftarrow\\beta M\+\(1\-\\beta\)G\_\{R\}
17:Approximately orthogonalize the buffered direction:

M~←NSOrth⁡\(M,TNS\)\\widetilde\{M\}\\leftarrow\\operatorname\{NSOrth\}\(M,T\_\{\\mathrm\{NS\}\}\)
18:Take a step and approximately project back to the Stiefel manifold:

Γ\+←Γ−η​M~,Γ←NSOrth⁡\(Γ\+,TNS\)\\Gamma\_\{\+\}\\leftarrow\\Gamma\-\\eta\\widetilde\{M\},\\qquad\\Gamma\\leftarrow\\operatorname\{NSOrth\}\(\\Gamma\_\{\+\},T\_\{\\mathrm\{NS\}\}\)
19:Convert

Γ\\Gammaback to the Kraus tensor

KK
20:endfor

21:if

DvalD\_\{\\mathrm\{val\}\}is providedthen

22:Evaluate validation log\-likelihood

ℓval\\ell\_\{\\mathrm\{val\}\}
23:if

ℓval\>ℓbest\\ell\_\{\\mathrm\{val\}\}\>\\ell\_\{\\mathrm\{best\}\}then

24:

ℓbest←ℓval\\ell\_\{\\mathrm\{best\}\}\\leftarrow\\ell\_\{\\mathrm\{val\}\}
25:

Kbest←KK\_\{\\mathrm\{best\}\}\\leftarrow K
26:endif

27:else

28:

Kbest←KK\_\{\\mathrm\{best\}\}\\leftarrow K
29:endif

30:Decay the stepsize:

η←α​η\\eta\\leftarrow\\alpha\\eta
31:endfor

32:return

KbestK\_\{\\mathrm\{best\}\}

Algorithm 2Newton–Schulz Approximate Orthogonalization1:Matrix

Z∈ℂp×nZ\\in\\mathbb\{C\}^\{p\\times n\}, number of iterations

TNST\_\{\\mathrm\{NS\}\}
2:Approximately orthogonalized matrix

ZZ
3:Set coefficients

a=3\.4445a=3\.4445,

b=−4\.7750b=\-4\.7750,

c=2\.0315c=2\.0315
4:Normalize

Z←Z/\(‖Z‖F\+ϵ\)Z\\leftarrow Z/\(\\\|Z\\\|\_\{F\}\+\\epsilon\)
5:if

p<np<nthen

6:

Z←Z†Z\\leftarrow Z^\{\\dagger\}
7:Mark that

ZZwas transposed

8:endif

9:for

t=1,…,TNSt=1,\\ldots,T\_\{\\mathrm\{NS\}\}do

10:

H←Z†​ZH\\leftarrow Z^\{\\dagger\}Z
11:

P←b​H\+c​H2P\\leftarrow bH\+cH^\{2\}
12:

Z←a​Z\+Z​PZ\\leftarrow aZ\+ZP
13:endfor

14:if

ZZwas transposedthen

15:

Z←Z†Z\\leftarrow Z^\{\\dagger\}
16:endif

17:return

ZZ

### 3\.2Mathematical Guarantees

We now state the assumptions and main mathematical guarantees for NS\-RIS\. AsBBis drawn i\.i\.d\. from the empirical distribution over the training set, we could, with a slight abuse of notation, write

F​\(Γ\)=𝔼B​\[L​\(Γ;B\)\],Γ∈St​\(p,n\),p=s​w​n\.F\(\\Gamma\)=\\mathbb\{E\}\_\{B\}\[L\(\\Gamma;B\)\],\\qquad\\Gamma\\in\\mathrm\{St\}\(p,n\),\\quad p=swn\.We use the real Frobenius inner product⟨Γ,Γ′⟩F=Re⁡tr⁡\(Γ†​Γ′\)\\langle\\Gamma,\\Gamma^\{\\prime\}\\rangle\_\{F\}=\\operatorname\{Re\}\\operatorname\{tr\}\(\\Gamma^\{\\dagger\}\\Gamma^\{\\prime\}\)\. The Riemannian gradient ofFFatΓ\\Gammais

grad⁡F​\(Γ\)=ΠΓ​\(∇F​\(Γ\)\),\\operatorname\{grad\}F\(\\Gamma\)=\\Pi\_\{\\Gamma\}\\bigl\(\\nabla F\(\\Gamma\)\\bigr\),where∇F​\(Γ\)\\nabla F\(\\Gamma\)denotes the Euclidean gradient andΠΓ​\(⋅\)\\Pi\_\{\\Gamma\}\(\\cdot\)is the orthogonal projection ontoTΓ​St​\(p,n\)T\_\{\\Gamma\}\\mathrm\{St\}\(p,n\)\(given in equation \([5](https://arxiv.org/html/2608.06554#S3.E5)\)\) under the real Frobenius inner product\.

For any matrixZZ, letPolar⁡\(Z\)\\operatorname\{Polar\}\(Z\)denote its thin polar factor\. That is, ifZ=U​Σ​V†Z=U\\Sigma V^\{\\dagger\}is a thin singular value decomposition, then

Polar⁡\(Z\)=U​V†and⟨Z,Polar⁡\(Z\)⟩F=‖Z‖1,\\operatorname\{Polar\}\(Z\)=UV^\{\\dagger\}\\quad\\text\{and\}\\quad\\langle Z,\\operatorname\{Polar\}\(Z\)\\rangle\_\{F\}=\\\|Z\\\|\_\{1\},where∥⋅∥1\\\|\\cdot\\\|\_\{1\}denotes the Schatten–11norm \(also known as the nuclear norm\)\. Throughout, we measure smoothness with respect to the Schatten–∞\\inftynorm,∥⋅∥∞\\\|\\cdot\\\|\_\{\\infty\}\(also known as the operator norm\), whose dual norm is the Schatten–11norm\. This operator–nuclear geometry is particularly natural for optimization algorithms based on polar factors or orthogonalized search directions\(Nesterov,,[2013](https://arxiv.org/html/2608.06554#bib.bib23); Jaggi,,[2013](https://arxiv.org/html/2608.06554#bib.bib14); Beck,,[2017](https://arxiv.org/html/2608.06554#bib.bib3)\), since

‖Z‖1=sup‖Y‖∞≤1⟨Z,Y⟩F=⟨Z,Polar⁡\(Z\)⟩F,\\\|Z\\\|\_\{1\}=\\sup\_\{\\\|Y\\\|\_\{\\infty\}\\leq 1\}\\langle Z,Y\\rangle\_\{F\}=\\langle Z,\\operatorname\{Polar\}\(Z\)\\rangle\_\{F\},where the supremum is attained atY=Polar⁡\(Z\)Y=\\operatorname\{Polar\}\(Z\)\.

###### Assumption 1\(Lipschitz smoothness\)\.

The objectiveFFis bounded below byF∗F^\{\*\}onSt​\(p,n\)\\mathrm\{St\}\(p,n\)\. Moreover,FFis continuously differentiable andLL\-Lipschitz smooth, i\.e\., for allΓ,Γ′∈St​\(p,n\)\\Gamma,\\Gamma^\{\\prime\}\\in\\mathrm\{St\}\(p,n\)there existsL\>0L\>0such that

‖∇F​\(Γ\)−∇F​\(Γ′\)‖1≤L​‖Γ−Γ′‖∞\\displaystyle\\\|\\nabla F\(\\Gamma\)\-\\nabla F\(\\Gamma^\{\\prime\}\)\\\|\_\{1\}\\;\\leq\\;L\\,\\\|\\Gamma\-\\Gamma^\{\\prime\}\\\|\_\{\\infty\}

Assumption[1](https://arxiv.org/html/2608.06554#Thmassumption1)is the standard smoothness condition used in nonconvex first\-order optimization\. We therefore measure stationarity using the Schatten\-11norm of the Riemannian gradient\. This criterion is stronger than the usual Frobenius\-norm criterion: since‖A‖F≤‖A‖1\\\|A\\\|\_\{F\}\\leq\\\|A\\\|\_\{1\}for every matrixAA, any point that isϵ\\epsilon\-stationary in nuclear norm is alsoϵ\\epsilon\-stationary in Frobenius norm\. In this paper, the convergence metric is defined as follows:

###### Definition 2\(ϵ\\epsilon\-stationary point\)\.

We callΓ∈ℝp×n\\Gamma\\in\\mathbb\{R\}^\{p\\times n\}anϵ\\epsilon\-stationary point \(in the nuclear norm\) if𝔼​\[‖grad⁡F​\(Γ\)‖1\]≤ϵ\\mathbb\{E\}\[\\\|\\operatorname\{grad\}F\(\\Gamma\)\\\|\_\{1\}\]\\leq\\epsilon\. Equivalently, we say an algorithm attainsϵ\\epsilon\-stationarity inTTsteps if

1T​∑t=1T𝔼​\[‖grad⁡F​\(Γt−1\)‖1\]≤ϵ\.\\displaystyle\\frac\{1\}\{T\}\\sum\_\{t=1\}^\{T\}\\mathbb\{E\}\[\\\|\\operatorname\{grad\}F\(\\Gamma\_\{t\-1\}\)\\\|\_\{1\}\]\\leq\\epsilon\.

Proposition[1](https://arxiv.org/html/2608.06554#Thmproposition1)explains the geometric role of the Newton–Schulz steps in NS\-RIS: the first step approximates the operator\-norm steepest descent direction through the polar factor of the momentum, while the second step acts as an efficient Stiefel retraction that restores feasibility after the update\. The proof is deferred to Appendix[A](https://arxiv.org/html/2608.06554#A1)\.

###### Proposition 1\.

Suppose Assumption[1](https://arxiv.org/html/2608.06554#Thmassumption1)holds\. LetMk\+1∈ℂp×nM\_\{k\+1\}\\in\\mathbb\{C\}^\{p\\times n\}be the momentum\-smoothed gradient used by NS\-RIS\. For a feasible first\-order perturbation atΓk\\Gamma\_\{k\}, represented by a curveγ\\gammaonSt​\(p,n\)\\mathrm\{St\}\(p,n\)with initial velocityΔ∈TΓk​St​\(p,n\)\\Delta\\in T\_\{\\Gamma\_\{k\}\}\\mathrm\{St\}\(p,n\), Lemma[1](https://arxiv.org/html/2608.06554#Thmlemma1)gives

F​\(γ​\(ηk\)\)=F​\(Γk\)\+ηk​⟨grad⁡F​\(Γk\),Δ⟩F\+O​\(ηk2\)\.F\(\\gamma\(\\eta\_\{k\}\)\)=F\(\\Gamma\_\{k\}\)\+\\eta\_\{k\}\\left\\langle\\operatorname\{grad\}F\(\\Gamma\_\{k\}\),\\Delta\\right\\rangle\_\{F\}\+O\(\\eta\_\{k\}^\{2\}\)\.Replacing the full Riemannian gradient by its momentum\-smoothed stochastic surrogateMk\+1M\_\{k\+1\}, NS\-RIS uses the corresponding operator\-norm linearized model\. The following ambient subproblem gives the steepest decrease direction for this model; feasibility of the iterate is restored by the retraction step below\. Consider

minΔ∈ℂp×n⟨Mk\+1,Δ⟩Fsubject to∥Δ∥∞≤1\.\\min\_\{\\Delta\\in\\mathbb\{C\}^\{p\\times n\}\}\\left\\langle M\_\{k\+1\},\\Delta\\right\\rangle\_\{F\}\\quad\\text\{subject to\}\\quad\\\|\\Delta\\\|\_\{\\infty\}\\leq 1\.IfMk\+1=U​Σ​V†M\_\{k\+1\}=U\\Sigma V^\{\\dagger\}is a thin singular value decomposition, then an optimizer is

Δ⋆=−U​V†=−Polar⁡\(Mk\+1\),\\Delta^\{\\star\}=\-UV^\{\\dagger\}=\-\\operatorname\{Polar\}\(M\_\{k\+1\}\),and the optimal value is−‖Mk\+1‖1\-\\\|M\_\{k\+1\}\\\|\_\{1\}\. Consequently, the first Newton–Schulz orthogonalization in Algorithm[1](https://arxiv.org/html/2608.06554#alg1),

M~k\+1=NSOrth⁡\(Mk\+1,TNS\),\\widetilde\{M\}\_\{k\+1\}=\\operatorname\{NSOrth\}\(M\_\{k\+1\},T\_\{\\mathrm\{NS\}\}\),computes a decomposition\-free approximation of the operator\-norm steepest descent direction\. The intermediate updateΓk\+12=Γk−ηk​M~k\+1\\Gamma\_\{k\+\\frac\{1\}\{2\}\}=\\Gamma\_\{k\}\-\\eta\_\{k\}\\widetilde\{M\}\_\{k\+1\}is then mapped back to the Stiefel manifold by the second Newton–Schulz orthogonalization\. In the exact polar projection case, this feasibility step solves

minQ∈St​\(p,n\)⁡‖Γk\+12−Q‖F2\\min\_\{Q\\in\\mathrm\{St\}\(p,n\)\}\\\|\\Gamma\_\{k\+\\frac\{1\}\{2\}\}\-Q\\\|\_\{F\}^\{2\}and has the form

R​\(Γk\+12\)=Γk\+12​\(Γk\+12†​Γk\+12\)−1/2,R\(\\Gamma\_\{k\+\\frac\{1\}\{2\}\}\)=\\Gamma\_\{k\+\\frac\{1\}\{2\}\}\\left\(\\Gamma\_\{k\+\\frac\{1\}\{2\}\}^\{\\dagger\}\\Gamma\_\{k\+\\frac\{1\}\{2\}\}\\right\)^\{\-1/2\},and satisfiesR​\(Γk\+12\)†​R​\(Γk\+12\)=𝕀nR\(\\Gamma\_\{k\+\\frac\{1\}\{2\}\}\)^\{\\dagger\}R\(\\Gamma\_\{k\+\\frac\{1\}\{2\}\}\)=\\mathbb\{I\}\_\{n\}wheneverΓk\+12\\Gamma\_\{k\+\\frac\{1\}\{2\}\}has full column rank\. The finite Newton–Schulz mapNSOrth⁡\(Γk\+12,TNS\)\\operatorname\{NSOrth\}\(\\Gamma\_\{k\+\\frac\{1\}\{2\}\},T\_\{\\mathrm\{NS\}\}\)is used as an efficient approximation of this projection/retraction\.

###### Assumption 2\(Stochastic projected gradients\)\.

At iterationkk, the sample indexiki\_\{k\}is drawn uniformly from the training set conditional on the historyℱk\\mathcal\{F\}\_\{k\}, and

gk=ΠΓk​\(∇ΓL​\(Γk;Bik\)\)g\_\{k\}=\\Pi\_\{\\Gamma\_\{k\}\}\\bigl\(\\nabla\_\{\\Gamma\}L\(\\Gamma\_\{k\};B\_\{i\_\{k\}\}\)\\bigr\)is an unbiased estimator of the Riemannian gradient:𝔼​\[gk∣ℱk\]=grad⁡F​\(Γk\)\\mathbb\{E\}\[g\_\{k\}\\mid\\mathcal\{F\}\_\{k\}\]=\\operatorname\{grad\}F\(\\Gamma\_\{k\}\)\. There is a constantσ2<∞\\sigma^\{2\}<\\inftysuch that

𝔼​\[‖gk−grad⁡F​\(Γk\)‖F2∣ℱk\]≤σ2\.\\mathbb\{E\}\\\!\\left\[\\\|g\_\{k\}\-\\operatorname\{grad\}F\(\\Gamma\_\{k\}\)\\\|\_\{F\}^\{2\}\\mid\\mathcal\{F\}\_\{k\}\\right\]\\leq\\sigma^\{2\}\.For a mini\-batch of sizebb, the right\-hand side is replaced byσ2/b\\sigma^\{2\}/b\.

Assumption[2](https://arxiv.org/html/2608.06554#Thmassumption2)is the standard unbiased bounded\-variance condition for stochastic gradient methods\. It follows when the mini\-batch samples are drawn independently from the empirical distribution and the per\-sample projected gradients have bounded second moment; on the compact Stiefel manifold this boundedness is automatic when the per\-sample gradients are continuous\. The factorσ2/b\\sigma^\{2\}/brecords the usual variance reduction from averagingbbi\.i\.d\. samples\. Assumption[3](https://arxiv.org/html/2608.06554#Thmassumption3)below only requires the finite Newton–Schulz orthogonalization to approximate the exact polar direction uniformly within a fixed tolerance; it does not require exact projection or an exact singular value decomposition\. This is mild in practice because Newton–Schulz iteration for approximating the polar factor is a classical and well\-studied approach, with convergence properties established in early work on iterative orthogonalization and polar decomposition and further refined by modern analyses of Newton\- and Halley\-type iterations for matrix polar decomposition and related matrix sign iterations\(Åke Björck and Bowie,,[1971](https://arxiv.org/html/2608.06554#bib.bib36); Higham,,[1986](https://arxiv.org/html/2608.06554#bib.bib11); Nakatsukasa et al\.,,[2010](https://arxiv.org/html/2608.06554#bib.bib21); Higham,,[2008](https://arxiv.org/html/2608.06554#bib.bib12)\)\.

###### Assumption 3\(Finite NS accuracy\)\.

Define the exact polar directionPk\+1=Polar⁡\(Mk\+1\)P\_\{k\+1\}=\\operatorname\{Polar\}\(M\_\{k\+1\}\)and the finite Newton–Schulz directionM~k\+1=NSOrth⁡\(Mk\+1,TNS\)\\widetilde\{M\}\_\{k\+1\}=\\operatorname\{NSOrth\}\(M\_\{k\+1\},T\_\{\\mathrm\{NS\}\}\)used in Algorithm[1](https://arxiv.org/html/2608.06554#alg1)\. For the chosen value ofTNST\_\{\\mathrm\{NS\}\}, there is constantεNS∈\[0,1\)\\varepsilon\_\{\\mathrm\{NS\}\}\\in\[0,1\)such that

‖M~k\+1−Pk\+1‖op≤εNS,\\\|\\widetilde\{M\}\_\{k\+1\}\-P\_\{k\+1\}\\\|\_\{\\mathrm\{op\}\}\\leq\\varepsilon\_\{\\mathrm\{NS\}\},for all iterations considered\.

Theorem[1](https://arxiv.org/html/2608.06554#Thmtheorem1)gives the finite\-time convergence guarantee for NS\-RIS\. The proof is deferred to Appendix[A](https://arxiv.org/html/2608.06554#A1)\.

###### Theorem 1\.

Suppose Assumptions[1](https://arxiv.org/html/2608.06554#Thmassumption1),[2](https://arxiv.org/html/2608.06554#Thmassumption2), and[3](https://arxiv.org/html/2608.06554#Thmassumption3)hold\. LetD=F​\(Γ0\)−F∗D=F\(\\Gamma\_\{0\}\)\-F^\{\*\}, let

G1≥supΓ∈St​\(p,n\)‖∇F​\(Γ\)‖1,LR=2​\(L\+G1\),G\_\{1\}\\geq\\sup\_\{\\Gamma\\in\\mathrm\{St\}\(p,n\)\}\\\|\\nabla F\(\\Gamma\)\\\|\_\{1\},\\qquad L\_\{R\}=2\(L\+G\_\{1\}\),and run Algorithm[1](https://arxiv.org/html/2608.06554#alg1)forKKinner iterations with mini\-batch sizebband a fixed numberTNST\_\{\\mathrm\{NS\}\}of Newton–Schulz iterations\. Let0≤β<10\\leq\\beta<1, letηk=η\>0\\eta\_\{k\}=\\eta\>0be constant, and define

Ek=Γk−\(Γk−1−η​M~k\)andρk=‖Ek‖∞\.E\_\{k\}=\\Gamma\_\{k\}\-\\bigl\(\\Gamma\_\{k\-1\}\-\\eta\\widetilde\{M\}\_\{k\}\\bigr\)\\quad\\text\{and\}\\quad\\rho\_\{k\}=\\\|E\_\{k\}\\\|\_\{\\infty\}\.For

ρ¯K=1K​∑k=1K𝔼​ρk,q¯K=1K​∑k=1K𝔼​ρk2,\\bar\{\\rho\}\_\{K\}=\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}\\mathbb\{E\}\\rho\_\{k\},\\qquad\\bar\{q\}\_\{K\}=\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\}\\mathbb\{E\}\\rho\_\{k\}^\{2\},and

𝒯K=1\(1−β\)​K​𝔼​‖M1−grad⁡F​\(Γ0\)‖1\+n​σb\+β​LR1−β​\(η​\(1\+εNS\)\+ρ¯K\)\.\\displaystyle\\mathcal\{T\}\_\{K\}=\\frac\{1\}\{\(1\-\\beta\)K\}\\mathbb\{E\}\\\|M\_\{1\}\-\\operatorname\{grad\}F\(\\Gamma\_\{0\}\)\\\|\_\{1\}\+\\frac\{\\sqrt\{n\}\\,\\sigma\}\{\\sqrt\{b\}\}\+\\frac\{\\beta L\_\{R\}\}\{1\-\\beta\}\\left\(\\eta\(1\+\\varepsilon\_\{\\mathrm\{NS\}\}\)\+\\bar\{\\rho\}\_\{K\}\\right\)\.the iterates satisfy

1K​∑k=0K−1𝔼​‖grad⁡F​\(Γk\)‖1≤\\displaystyle\\frac\{1\}\{K\}\\sum\_\{k=0\}^\{K\-1\}\\mathbb\{E\}\\\|\\operatorname\{grad\}F\(\\Gamma\_\{k\}\)\\\|\_\{1\}\\leqDη​\(1−εNS\)​K\+LR​\(1\+εNS\)2​η1−εNS\+2​G1​ρ¯K\+LR​q¯Kη​\(1−εNS\)\+2​𝒯K1−εNS\.\\displaystyle\\frac\{D\}\{\\eta\(1\-\\varepsilon\_\{\\mathrm\{NS\}\}\)K\}\+\\frac\{L\_\{R\}\(1\+\\varepsilon\_\{\\mathrm\{NS\}\}\)^\{2\}\\eta\}\{1\-\\varepsilon\_\{\\mathrm\{NS\}\}\}\+\\frac\{2G\_\{1\}\\bar\{\\rho\}\_\{K\}\+L\_\{R\}\\bar\{q\}\_\{K\}\}\{\\eta\(1\-\\varepsilon\_\{\\mathrm\{NS\}\}\)\}\+\\frac\{2\\mathcal\{T\}\_\{K\}\}\{1\-\\varepsilon\_\{\\mathrm\{NS\}\}\}\.

## 4Numerical Analysis

This section describes the experimental setup and evaluation metric in Subsection[4\.1](https://arxiv.org/html/2608.06554#S4.SS1), then reports performance on the synthetic HMM benchmark in Subsection[4\.2](https://arxiv.org/html/2608.06554#S4.SS2)and the synthetic HQMM benchmark in Subsection[4\.3](https://arxiv.org/html/2608.06554#S4.SS3)\.

### 4\.1Setup and Evaluation Metric

We use the experimental setting and evaluation convention as inAdhikary et al\., \([2020](https://arxiv.org/html/2608.06554#bib.bib1)\), which is a standard benchmark protocol for HQMM learning algorithms\. For an HQMM with latent dimensionnn, output alphabet sizess, andwwKraus operators per output, we write the model as an\(n,s,w\)\(n,s,w\)\-HQMM\. Equivalently, after stacking all Kraus operators vertically, the optimization variable is a complex Stiefel matrixΓ∈ℂn​s​w×n\\Gamma\\in\\mathbb\{C\}^\{nsw\\times n\}whose columns satisfyΓ†​Γ=𝕀n\\Gamma^\{\\dagger\}\\Gamma=\\mathbb\{I\}\_\{n\}\.

Unless otherwise stated, the latent density matrix is initialized as a random Hermitian positive semidefinite matrix with unit trace, and the stacked Kraus matrix is initialized as a random orthonormal matrix\. Training minimizes the negative log\-likelihood of the observed sequences under the HQMM\. During optimization, validation performance is monitored after each epoch, and the parameter tensor with the largest validation log\-likelihood is retained for final evaluation\. This selection rule is equivalent to choosing the model that assigns the highest probability to the validation set, while leaving the reported test metric independent of the particular sequence lengths used during training\.

We compare against three standard learning procedures\. The Givens Search \(GS\) method of[Srinivasan et al\., 2018b](https://arxiv.org/html/2608.06554#bib.bib30)maintains feasibility by applying local Givens rotations, or unitary transformations, to the stacked Kraus matrix and accepting likelihood\-improving updates\. The constrained optimization on the Stiefel manifold \(COSM\) method ofAdhikary et al\., \([2020](https://arxiv.org/html/2608.06554#bib.bib1)\)instead performs gradient\-based updates directly on the Stiefel manifold using a retraction that we propose\. In the experiments ofAdhikary et al\., \([2020](https://arxiv.org/html/2608.06554#bib.bib1)\), COSM converged to better optima faster than GS and scaled to larger HQMMs that were too slow to train with GS\. The Expectation\-Maximization \(EM\) algorithm is the standard maximum\-likelihood procedure for HMMs, alternating between inference of latent\-state responsibilities and parameter re\-estimation; we use EM\-trained HMMs as the classical baseline\. On the synthetic HMM benchmark,Adhikary et al\., \([2020](https://arxiv.org/html/2608.06554#bib.bib1)\)report that small HQMMs can outperform small HMMs, although this advantage does not hold for the66\-state HMM baseline, so EM remains a strong reference method when the data are generated by an HMM\.

In this section, we evaluate the synthetic HMM and HQMM benchmarks\. For the HMM benchmark, we follow the synthetic data setting of[Srinivasan et al\., 2018b](https://arxiv.org/html/2608.06554#bib.bib30)andAdhikary et al\., \([2020](https://arxiv.org/html/2608.06554#bib.bib1)\): the data are generated by an HMM with66hidden states and66possible outputs, and we use the same2020training and1010validation sequences of length30003000as they did\. Each long sequence is split into300300shorter sequences, with a burn\-in of100100, before training HQMMs for6060epochs and selecting the model with the highest validation description accuracy\. For the HQMM benchmark, we use the synthetic HQMM data in those two papers, generated by a22\-hidden\-state,66\-output HQMM\. This benchmark uses the same2020training and1010validation sequences of length30003000, the same split into300300shorter sequences, and the same burn\-in of100100for training; evaluation is performed on1010test sequences of length30003000using a burn\-in of10001000\. This shortened training protocol reduces run time without changing the amount of training data processed\.

For all these experiments, we report description accuracy, a scaled log\-likelihood metric used in prior HQMM work\([Srinivasan et al\., 2018b,](https://arxiv.org/html/2608.06554#bib.bib30); Adhikary et al\.,,[2020](https://arxiv.org/html/2608.06554#bib.bib1)\)\. Given a test sequenceYYof lengthℓ\\ell, an alphabet of sizess, and a trained model𝔻\\mathbb\{D\}, the description accuracy is

D​A=f​\(1\+logs⁡P​\(Y∣𝔻\)ℓ\),DA=f\\left\(1\+\\frac\{\\log\_\{s\}P\(Y\\mid\\mathbb\{D\}\)\}\{\\ell\}\\right\),where

f​\(x\)=\{tanh⁡\(x/8\),x≤0,x,x\>0\.f\(x\)=\\begin\{cases\}\\tanh\(x/8\),&x\\leq 0,\\\\ x,&x\>0\.\\end\{cases\}The normalization byℓ\\ellmakes the score comparable across sequences of different lengths\. A value ofD​A=1DA=1corresponds to assigning probability one to the observed sequence, whileD​A\>0DA\>0indicates performance better than the uniform random baseline over the output alphabet\. Thus, higher values ofD​ADAindicate better predictive performance\. When multiple test sequences are evaluated, we report the mean description accuracy and use the standard deviation across test sequences as the error bar\. For labeled sequence\-classification experiments, we follow the same likelihood\-based decision rule as in HQMM benchmarks: one model is trained for each class, and the predicted label is the class whose model assigns the highest likelihood to the test sequence; the reported metric is average classification accuracy\.

### 4\.2Synthetic HMM Test Performance

Figure[1](https://arxiv.org/html/2608.06554#S4.F1)compares HQMM learning methods under different hidden dimensions and Kraus ranks\. In the left panel, we fix the Kraus rank tow=1w=1and vary the hidden dimensionnn\. Under this setting, NS\-RIS consistently achieves the highest test description accuracy among all HQMM training methods\. The advantage becomes more pronounced as the hidden dimension increases, with NS\-RIS substantially outperforming both GS and COSM for largernn, while also exceeding the corresponding classical EM baseline\. In contrast, GS and COSM fail to surpass EM whenn=6n=6\. In the right panel, we fix the hidden dimension atn=6n=6and vary the Kraus rankwwto evaluate the benefit of increasing the number of Kraus operators\. Since EM does not use Kraus operators, its performance depends only on\(n,s\)\(n,s\)and is therefore independent ofww\. NS\-RIS maintains a clear performance advantage across all tested Kraus ranks, demonstrating that the proposed update rule can effectively exploit the richer HQMM parameterization\. GS results are omitted forw\>1w\>1, because the method becomes computationally prohibitive in these cases as revealed inAdhikary et al\., \([2020](https://arxiv.org/html/2608.06554#bib.bib1)\)\.

![Refer to caption](https://arxiv.org/html/2608.06554v1/pic/hmm_combined.png)Figure 1:Test description accuracy for HQMM learning methods under varying hidden dimensions and Kraus ranks\. Left: Performance with Kraus rank fixed atw=1w=1while varying the hidden dimensionnn\. We report the mean test description accuracy, with error bars indicating the standard deviation across test sequences\. Right: Performance for models with hidden dimensionn=6n=6while varying the Kraus rankww\. The EM baseline is shown as a horizontal dashed line because its performance is independent ofww\.In Figure[2](https://arxiv.org/html/2608.06554#S4.F2), we compare the running times of the HQMM learning methods under varying hidden dimensions and Kraus ranks\. In the left panel, we fix the Kraus rank atw=1w=1and vary the hidden dimensionnn\. NS\-RIS and COSM exhibit nearly identical running times across all tested hidden dimensions, remaining relatively stable asnnincreases\. In contrast, GS incurs substantially higher computational cost, with runtime increasing steadily as the hidden dimension grows\. This highlights the significantly better scalability of NS\-RIS and COSM with respect to the latent dimension\. In the right panel, we fix the hidden dimension atn=6n=6and vary the Kraus rankww\. Both NS\-RIS and COSM experience increased computational cost as the Kraus rank grows, reflecting the larger HQMM parameterization\. However, NS\-RIS remains competitive with COSM across all tested Kraus ranks and is slightly faster for larger values ofww\. Combined with the accuracy results in Figure[1](https://arxiv.org/html/2608.06554#S4.F1), these findings demonstrate that NS\-RIS achieves stronger predictive performance without incurring additional computational overhead relative to existing scalable HQMM training methods\.

![Refer to caption](https://arxiv.org/html/2608.06554v1/pic/hmm_runtime.png)Figure 2:Running time comparison for HQMM learning methods under varying hidden dimensions and Kraus ranks\. Left: Running time with Kraus rank fixed atw=1w=1while varying the hidden dimensionnn\. NS\-RIS and COSM maintain relatively stable runtimes across hidden dimensions, whereas GS becomes substantially more expensive asnnincreases\. Right: Running time for models with hidden dimensionn=6n=6while varying the Kraus rankww\. Both NS\-RIS and COSM require additional computation for larger Kraus ranks, but NS\-RIS remains competitive with COSM across all tested settings\.Across the evaluated synthetic benchmark settings, NS\-RIS consistently outperforms COSM in terms of test performance\. To quantify the relative gain, we compute the percentage improvement using

Improvement\(%\)=Metric∗NS\-RIS−Metric∗COSMMetricCOSM×100\.\\mathrm\{Improvement\}\(\\%\)=\\frac\{\\mathrm\{Metric\}\*\{\\text\{NS\-RIS\}\}\-\\mathrm\{Metric\}\*\{\\text\{COSM\}\}\}\{\\mathrm\{Metric\}\_\{\\text\{COSM\}\}\}\\times 100\.Averaged across all experimental configurations, NS\-RIS achieves approximately38\.5%38\.5\\%relative improvement over COSM\. The largest improvement is observed forscenario\_id=4\\texttt\{scenario\\\_id\}=4, where the averaged test metric increases from0\.17360\.1736for COSM to0\.26150\.2615for NS\-RIS\. Substituting these values into the relative\-improvement formula gives0\.2615−0\.17360\.1736×100=50\.6%\\frac\{0\.2615\-0\.1736\}\{0\.1736\}\\times 100=50\.6\\%, demonstrating that NS\-RIS attains over50%50\\%relative improvement in the best\-performing configuration\. These results indicate that the proposed NS\-RIS optimization framework provides substantial and consistent performance gains over existing geometric optimization approaches for synthetic HMM learning tasks\.

It is worth noting that only approximately10%10\\%of the experimental settings involved explicit hyperparameter tuning\. Specifically, hyperparameter tuning was conducted only for the hidden\-state configuration withn=4n=4, while for all remaining experiments we directly adopted the best\-performing hyperparameter settings reported for COSM inAdhikary et al\., \([2020](https://arxiv.org/html/2608.06554#bib.bib1)\)\. Despite this minimal tuning effort, NS\-RIS consistently achieved superior predictive performance and competitive runtime efficiency across the evaluated benchmarks\. This result highlights the robustness and strong generalization capability of the proposed NS\-RIS optimization framework, suggesting that it is substantially less sensitive to hyperparameter selection while still converging to high\-quality solutions\. The same phenomenon was observed in the subsequent HQMM benchmark experiments\.

### 4\.3Synthetic HQMM Test Performance

The synthetic HQMM benchmark follows the experimental setting introduced byAdhikary et al\., \([2020](https://arxiv.org/html/2608.06554#bib.bib1)\), where the data are generated from a synthetic HQMM inspired by the Stern–Gerlach experiment in quantum mechanics\. The benchmark uses a configuration with hidden dimensionn=2n=2, output alphabet sizes=6s=6, and Kraus rankw=1w=1, which was specifically designed to highlight the expressive advantages of HQMMs over classical HMMs\. In the original study, the authors demonstrated that a significantly larger classical HMM is required to match the representational capacity of a small HQMM, emphasizing the richer expressiveness of quantum\-inspired latent\-state models\. Figure[3](https://arxiv.org/html/2608.06554#S4.F3)presents the averaged training, validation, and test metrics together with standard\-deviation error bars across repeated runs, as well as the running\-time comparison between NS\-RIS, GS, COSM, and EM\. The results show that NS\-RIS consistently achieves the strongest predictive performance across all evaluation metrics while maintaining competitive computational efficiency\. In particular, NS\-RIS improves the test metric over COSM by approximately18\.9%18\.9\\%, computed as0\.1536−0\.12920\.1292×100\\frac\{0\.1536\-0\.1292\}\{0\.1292\}\\times 100, while simultaneously reducing runtime by approximately12\.0%12\.0\\%relative to COSM, computed as643\.23−566\.19643\.23×100\\frac\{643\.23\-566\.19\}\{643\.23\}\\times 100\. These results demonstrate that the proposed NS\-RIS can improve both solution quality and practical scalability for HQMM learning\. It is worth mentioning that no additional hyperparameter tuning was conducted in our experiments\. All results were generated using the same hyperparameter configuration reported as optimal for the COSM method inAdhikary et al\., \([2020](https://arxiv.org/html/2608.06554#bib.bib1)\), and this identical configuration was applied uniformly across NS\-RIS, GS, COSM, and EM to ensure a fair and controlled comparison between optimization methods\.

![Refer to caption](https://arxiv.org/html/2608.06554v1/pic/hqmm_one_scenario_comparison.png)Figure 3:Comparison of NS\-RIS, GS, COSM, and EM on the synthetic HQMM benchmark with configuration\(n=2,s=6,w=1\)\(n=2,s=6,w=1\)\. The top row reports the averaged training, validation, and test metrics with standard\-deviation error bars across repeated runs, while the bottom row compares elapsed running times on a logarithmic scale\. NS\-RIS consistently achieves the best predictive performance across all evaluation metrics while also providing competitive computational efficiency\.

## 5Empirical Analysis

This section evaluates NS\-RIS on the real Splice benchmark\. We first describe the dataset and preprocessing in Subsection[5\.1](https://arxiv.org/html/2608.06554#S5.SS1), then report the classification and runtime results in Subsection[5\.2](https://arxiv.org/html/2608.06554#S5.SS2), and finally discuss the biological explanation for using higher\-dimensional latent spaces in Subsection[5\.3](https://arxiv.org/html/2608.06554#S5.SS3)\.

### 5\.1The Splice Dataset

The splice dataset is a widely used benchmark for sequence classification and splice\-junction prediction in computational biology\(Towell and Shavlik,,[1991](https://arxiv.org/html/2608.06554#bib.bib31); Dheeru and Karra Taniskidou,,[2017](https://arxiv.org/html/2608.06554#bib.bib8)\)that is publicly available at UCI machine learning repository\. The dataset consists of DNA sequences of fixed length6060, where each position corresponds to one of the four nucleobases: Adenine \(A\), Cytosine \(C\), Guanine \(G\), or Thymine \(T\)\. Biologically, DNA sequences contain protein\-coding regions known as exons, interspersed with non\-coding regions called introns\. Correct identification of exon–intron boundaries is an important problem in gene prediction and genomic analysis\. The classification task is to determine whether a sequence contains an exon–intron \(EI\) splice site, an intron–exon \(IE\) splice site, or neither \(N\)\. The dataset contains762762EI examples,765765IE examples, and16481648negative examples\. Although the original dataset includes several ambiguous nucleotide symbols beyondA,C,G,T\{A,C,G,T\}, we remove sequences containing ambiguous characters during preprocessing in order to ensure a consistent discrete input alphabet for training and evaluation\.

### 5\.2Splice Classification Results

Figure[4](https://arxiv.org/html/2608.06554#S5.F4)compares mean classification error for COSM and NS\-RIS over latent dimensionsn∈\{2,4,6,8\}n\\in\\\{2,4,6,8\\\}and Kraus ranksw∈\{1,2,4,6\}w\\in\\\{1,2,4,6\\\}\. The results show that increasing the latent dimension is beneficial for NS\-RIS: forn=6n=6andn=8n=8, NS\-RIS consistently achieves lower error than the corresponding COSM configurations and falls well below the EM baseline, with the best mean error reaching approximately0\.3400\.340forn=6n=6and0\.3300\.330forn=8n=8\. More specifically, atn=6n=6, the best COSM error is0\.4220\.422, whereas NS\-RIS reaches0\.3400\.340, a relative reduction of about19\.4%19\.4\\%; this is also substantially below the EM baseline error0\.4280\.428\. Atn=8n=8, the best COSM error is0\.4020\.402, while NS\-RIS achieves0\.3300\.330, a relative reduction of about17\.9%17\.9\\%and an improvement over the EM baseline error0\.4210\.421\. For smaller latent dimensions \(n=2n=2andn=4n=4\), the performance is more mixed, indicating that the advantage of NS\-RIS becomes most visible when the HQMM has enough latent capacity to model the splice\-junction structure\.

![Refer to caption](https://arxiv.org/html/2608.06554v1/pic/splice_combined_mean_error_bars.png)Figure 4:Mean error comparison of COSM and NS\-RIS on the Splice dataset for latent state sizesn∈\{2,4,6,8\}n\\in\\\{2,4,6,8\\\}\. Within each panel, bars report the mean error for window sizesw∈\{1,2,4,6\}w\\in\\\{1,2,4,6\\\}, with darker shades indicating largerww\. The dashed horizontal line shows the EM baseline for the same latent state size:0\.4580\.458forn=2n=2,0\.4560\.456forn=4n=4,0\.4280\.428forn=6n=6, and0\.4210\.421forn=8n=8\.Figure[5](https://arxiv.org/html/2608.06554#S5.F5)gives a class\-wise comparison between COSM and NS\-RIS using the EI, IE, and negative\-class error metrics\. The nested pies show that NS\-RIS improves all three class\-specific errors in the larger\-latent\-state settings\(n,w\)=\(6,4\)\(n,w\)=\(6,4\)and\(8,4\)\(8,4\), with particularly large reductions in the negative\-class error\. For\(n,w\)=\(6,4\)\(n,w\)=\(6,4\), NS\-RIS reduces the EI error from0\.4221600\.422160to0\.3490460\.349046\(17\.32%17\.32\\%lower\), the IE error from0\.3571890\.357189to0\.3284320\.328432\(8\.05%8\.05\\%lower\), and the negative\-class error from0\.5446040\.544604to0\.4094220\.409422\(24\.82%24\.82\\%lower\)\. For\(n,w\)=\(8,4\)\(n,w\)=\(8,4\), NS\-RIS reduces the EI error from0\.3897340\.389734to0\.3401660\.340166\(12\.72%12\.72\\%lower\), the IE error from0\.3467320\.346732to0\.3447710\.344771\(0\.57%0\.57\\%lower\), and the negative\-class error from0\.5307870\.530787to0\.3930540\.393054\(25\.95%25\.95\\%lower\)\. In the smaller settings, NS\-RIS is not uniformly better across all classes: for\(2,4\)\(2,4\)it improves IE error from0\.3888890\.388889to0\.3114380\.311438but increases EI and negative\-class errors, while for\(4,4\)\(4,4\)all three errors are higher than COSM\. This class\-wise view reinforces that the strongest advantage appears in the larger latent\-state regimes, where NS\-RIS improves not only the average error but also all three biologically meaningful classification categories\.

![Refer to caption](https://arxiv.org/html/2608.06554v1/pic/cosm_vs_nsris_three_metrics.png)Figure 5:Nested pie comparison of COSM and NS\-RIS errors across three metrics\.Figure[6](https://arxiv.org/html/2608.06554#S5.F6)compares the running time of NS\-RIS and COSM across the tested Splice configurations\. The upper subfigure shows that COSM generally becomes slower as the Kraus rankwwincreases, although its ordering is not perfectly monotone across all latent dimensions\. The lower subfigure shows a clearer monotone pattern for NS\-RIS: the running time grows first with the Kraus rankwwand then with the combined scalen​wnw, which more fully reflects the dimension of the stacked Stiefel parameter\. The runtime advantage of NS\-RIS is especially clear for several larger Kraus\-rank settings where COSM requires substantially more time\. Thus the improved classification accuracy is not obtained by paying a larger computational cost; rather, the Newton–Schulz retraction strategy improves predictive performance while preserving practical scalability\.

![Refer to caption](https://arxiv.org/html/2608.06554v1/pic/splice_runtime_nsris_cosm_bars.png)Figure 6:Running\-time comparison of NS\-RIS and COSM on the Splice dataset across latent dimensions and Kraus ranks\.
### 5\.3Scientific Validation

Recent advances in quantum computing and quantum machine learning have stimulated growing interest in applying quantum\-inspired methods to genomic and DNA sequence analysis\. Several studies have explored quantum algorithms for biological sequence comparison, DNA sequence alignment, genome assembly, and genomic pattern recognition\(Outeiral et al\.,,[2021](https://arxiv.org/html/2608.06554#bib.bib25); Kösoglu\-Kind et al\.,,[2023](https://arxiv.org/html/2608.06554#bib.bib16); Varsamis et al\.,,[2023](https://arxiv.org/html/2608.06554#bib.bib32)\)\. Broader surveys of quantum computing in biology likewise identify genetics, molecular biology, drug design, and bioinformatics as application areas where quantum representations may become useful\(Ghamsari,,[2025](https://arxiv.org/html/2608.06554#bib.bib10)\)\. In particular, quantum\-based approaches have been proposed to model high\-dimensional sequence dependencies and complex biological correlations that are difficult to capture using classical probabilistic models alone\(Nalecz\-Charkiewicz et al\.,,[2024](https://arxiv.org/html/2608.06554#bib.bib22)\)\. Recent work has also demonstrated the feasibility of encoding and analyzing genomic sequences on real quantum hardware, highlighting the potential of quantum representations for large\-scale genomic inference tasks\. These developments suggest that richer latent\-state frameworks, such as HQMMs, may provide advantages over classical HMMs in modeling splice\-junction sequences, especially when long\-range dependencies and higher\-order contextual interactions are present in biological sequence data\.

Although the observed alphabet in the splice dataset consists of only four nucleotidesA,C,G,T\{A,C,G,T\}, there is strong evidence from computational biology and sequence modeling that the latent state space should generally be substantially larger than four\. In HMMs and related state\-space approaches for genomic sequences, the latent variables are not intended to represent the nucleotides themselves, but rather higher\-order biological contexts such as splice donor/acceptor motifs, exon and intron regions, insertion/deletion dynamics, codon structure, and long\-range dependencies\. Classical profile\-HMM formulations for biological sequence analysis therefore employ many more hidden states than observable symbols, including match, insertion, and deletion states at different sequence positions\(Krogh et al\.,,[1994](https://arxiv.org/html/2608.06554#bib.bib17); Eddy,,[1998](https://arxiv.org/html/2608.06554#bib.bib9)\)\. In splice\-site prediction specifically, different latent states are used to model biologically distinct regimes surrounding exon–intron boundaries, even when the emissions remain limited to four nucleobases\(Burge and Karlin,,[1997](https://arxiv.org/html/2608.06554#bib.bib5)\)\. Consequently, restricting the latent dimensionality to four would force the model to conflate multiple heterogeneous biological mechanisms into the same representation, reducing its ability to capture contextual sequence structure\.

More recent representation\-learning approaches for DNA sequence modeling further support the use of latent spaces whose dimensionality exceeds the cardinality of the nucleotide alphabet\. Deep latent\-variable and sequence\-embedding models trained on splice\-site datasets consistently learn distributed representations that encode motif composition, positional dependencies, and regulatory patterns that cannot be represented using only four discrete latent categories\(Agarwal et al\.,,[2019](https://arxiv.org/html/2608.06554#bib.bib2); Poddar et al\.,,[2023](https://arxiv.org/html/2608.06554#bib.bib26)\)\. From an information\-theoretic perspective, the observable alphabet size only constrains the emission space, whereas the latent state dimension reflects the complexity of the underlying generative process\. Since splice recognition depends on combinatorial sequence motifs and contextual interactions extending across many nucleotide positions, using a latent dimension larger than four is both biologically motivated and empirically supported in prior literature\.

## 6Conclusion

This paper introduces NS\-RIS for scalable learning HQMM parameters on the complex Stiefel manifold\. NS\-RIS enables HQMM training to be both geometrically faithful and computationally efficient\. Specifically, the first Newton–Schulz iteration approximates the operator\-norm steepest descent direction via the polar factor of the momentum, while the second Newton–Schulz iteration restores the Stiefel\-manifold feasibility required by the trace\-preserving Kraus constraint\.

We establish the theoretical convergence properties of NS\-RIS under standard assumptions on smoothness, stochastic gradients, and finite Newton–Schulz accuracy\. In particular, we derive an explicit finite\-time stationarity bound that separates the contributions of the telescoping descent term, smoothness error, stochastic momentum\-tracking error, and feasibility residuals, thereby clarifying how each component influences convergence\. The analysis further shows that sufficiently accurate Newton–Schulz retraction errors contribute only lower\-order terms, implying that the decomposition\-free implementation retains the convergence behavior of stochastic Riemannian gradient descent up to the inherent stochastic noise floor\.

Extensive experiments on both synthetic and real\-world sequence datasets demonstrate the effectiveness of NS\-RIS\. On synthetic HMM and HQMM benchmarks, NS\-RIS consistently achieves stronger predictive performance than existing HQMM training methods while remaining computationally competitive\. On the Splice benchmark, NS\-RIS is particularly effective when the latent dimension is sufficiently large to capture the biological complexity of splice\-junction sequences, reducing both the overall classification error and class\-specific errors compared with the current state\-of\-the\-art method in higher\-dimensional settings\.

Importantly, NS\-RIS provides the first evidence in these benchmarks that an HQMM can significantly outperform the HMM baseline on data that are not generated by a quantum model\. This is a remarkable step beyond the theoretical fact that HQMMs generalize HMMs: it shows that, when equipped with scalable Stiefel\-manifold inference, HQMMs can deliver practical gains on ordinary sequence data\. These results therefore position HQMMs as viable replacements or extensions of HMMs across the broad scientific sequence\-modeling settings where HMMs have long been standard tools and richer latent dynamics are needed\.

## Appendix AProofs

This appendix gives the proofs of the mathematical statements in Subsection[3\.2](https://arxiv.org/html/2608.06554#S3.SS2)\. Subsection[A\.1](https://arxiv.org/html/2608.06554#A1.SS1)proves Proposition[1](https://arxiv.org/html/2608.06554#Thmproposition1), and Subsection[A\.2](https://arxiv.org/html/2608.06554#A1.SS2)proves Theorem[1](https://arxiv.org/html/2608.06554#Thmtheorem1)together with the auxiliary lemmas used in the convergence analysis\.

### A\.1Proof of Proposition[1](https://arxiv.org/html/2608.06554#Thmproposition1)

Lemma[1](https://arxiv.org/html/2608.06554#Thmlemma1)below gives the first\-order expansion on the Stiefel manifold\.

###### Lemma 1\.

LetFFbe continuously differentiable in a neighborhood ofSt​\(p,n\)\\mathrm\{St\}\(p,n\), and letΓ∈St​\(p,n\)\\Gamma\\in\\mathrm\{St\}\(p,n\)\. Ifγ:\(−ϵ,ϵ\)→St​\(p,n\)\\gamma:\(\-\\epsilon,\\epsilon\)\\to\\mathrm\{St\}\(p,n\)is a differentiable curve satisfyingγ​\(0\)=Γ\\gamma\(0\)=\\Gammaandγ′​\(0\)=ξ∈TΓ​St​\(p,n\)\\gamma^\{\\prime\}\(0\)=\\xi\\in T\_\{\\Gamma\}\\mathrm\{St\}\(p,n\), then

F​\(γ​\(t\)\)=F​\(Γ\)\+t​⟨grad⁡F​\(Γ\),ξ⟩F\+o​\(t\)\.F\(\\gamma\(t\)\)=F\(\\Gamma\)\+t\\left\\langle\\operatorname\{grad\}F\(\\Gamma\),\\xi\\right\\rangle\_\{F\}\+o\(t\)\.If the gradient ofFFis locally Lipschitz, the remainder isO​\(t2\)O\(t^\{2\}\); in particular, this quadratic remainder holds under Assumption[1](https://arxiv.org/html/2608.06554#Thmassumption1)\. In particular, for a straight perturbationΓ\+t​Δ\\Gamma\+t\\Delta, the same first\-order expression withgrad⁡F​\(Γ\)\\operatorname\{grad\}F\(\\Gamma\)is valid whenΔ∈TΓ​St​\(p,n\)\\Delta\\in T\_\{\\Gamma\}\\mathrm\{St\}\(p,n\)\. For a general ambient perturbationΔ∈ℂp×n\\Delta\\in\\mathbb\{C\}^\{p\\times n\}, the linear term is instead⟨∇F​\(Γ\),Δ⟩F\\langle\\nabla F\(\\Gamma\),\\Delta\\rangle\_\{F\}\.

###### Proof\.

By the chain rule and the definition of the Euclidean gradient,

dd​t​F​\(γ​\(t\)\)\|t=0=⟨∇F​\(Γ\),ξ⟩F\.\\frac\{d\}\{dt\}F\(\\gamma\(t\)\)\\bigg\|\_\{t=0\}=\\left\\langle\\nabla F\(\\Gamma\),\\xi\\right\\rangle\_\{F\}\.Sincegrad⁡F​\(Γ\)=ΠΓ​\(∇F​\(Γ\)\)\\operatorname\{grad\}F\(\\Gamma\)=\\Pi\_\{\\Gamma\}\(\\nabla F\(\\Gamma\)\)is the orthogonal projection of∇F​\(Γ\)\\nabla F\(\\Gamma\)ontoTΓ​St​\(p,n\)T\_\{\\Gamma\}\\mathrm\{St\}\(p,n\), the residual∇F​\(Γ\)−grad⁡F​\(Γ\)\\nabla F\(\\Gamma\)\-\\operatorname\{grad\}F\(\\Gamma\)is orthogonal to every tangent vector\. Becauseξ∈TΓ​St​\(p,n\)\\xi\\in T\_\{\\Gamma\}\\mathrm\{St\}\(p,n\),

⟨∇F​\(Γ\),ξ⟩F=⟨grad⁡F​\(Γ\),ξ⟩F\.\\left\\langle\\nabla F\(\\Gamma\),\\xi\\right\\rangle\_\{F\}=\\left\\langle\\operatorname\{grad\}F\(\\Gamma\),\\xi\\right\\rangle\_\{F\}\.This proves the stated first\-order expansion\. The first\-ordero​\(t\)o\(t\)term only uses differentiability\. If the gradient is locally Lipschitz, the usual Taylor estimate gives a quadratic remainder; this condition is supplied by Assumption[1](https://arxiv.org/html/2608.06554#Thmassumption1)in the operator–nuclear geometry used here\. The final claim follows by applying the same calculation to the straight curvet↦Γ\+t​Δt\\mapsto\\Gamma\+t\\Delta; whenΔ\\Deltais not tangent, the normal component ofΔ\\Deltaneed not be orthogonal to∇F​\(Γ\)\\nabla F\(\\Gamma\), so the Euclidean gradient gives the correct ambient first variation\. ∎

###### Proof of Proposition[1](https://arxiv.org/html/2608.06554#Thmproposition1)\.

By spectral–nuclear duality, for every feasibleΔ\\Delta,

⟨Mk\+1,Δ⟩F≥−‖Mk\+1‖1​‖Δ‖∞≥−‖Mk\+1‖1\.\\left\\langle M\_\{k\+1\},\\Delta\\right\\rangle\_\{F\}\\geq\-\\\|M\_\{k\+1\}\\\|\_\{1\}\\\|\\Delta\\\|\_\{\\infty\}\\geq\-\\\|M\_\{k\+1\}\\\|\_\{1\}\.TakingΔ⋆=−U​V†\\Delta^\{\\star\}=\-UV^\{\\dagger\}gives‖Δ⋆‖∞=1\\\|\\Delta^\{\\star\}\\\|\_\{\\infty\}=1and

⟨Mk\+1,Δ⋆⟩F=−⟨U​Σ​V†,U​V†⟩F=−tr⁡\(Σ\)=−‖Mk\+1‖1\.\\left\\langle M\_\{k\+1\},\\Delta^\{\\star\}\\right\\rangle\_\{F\}=\-\\left\\langle U\\Sigma V^\{\\dagger\},UV^\{\\dagger\}\\right\\rangle\_\{F\}=\-\\operatorname\{tr\}\(\\Sigma\)=\-\\\|M\_\{k\+1\}\\\|\_\{1\}\.Thus−Polar⁡\(Mk\+1\)\-\\operatorname\{Polar\}\(M\_\{k\+1\}\)solves the constrained linearized decrease problem\. For the projection statement, letY=Γk\+12Y=\\Gamma\_\{k\+\\frac\{1\}\{2\}\}andR​\(Y\)=Y​\(Y†​Y\)−1/2R\(Y\)=Y\(Y^\{\\dagger\}Y\)^\{\-1/2\}\. ForQ∈St​\(p,n\)Q\\in\\mathrm\{St\}\(p,n\),

‖Y−Q‖F2=tr⁡\(Y†​Y\)\+n−2​Re⁡tr⁡\(Q†​Y\)\.\\\|Y\-Q\\\|\_\{F\}^\{2\}=\\operatorname\{tr\}\(Y^\{\\dagger\}Y\)\+n\-2\\operatorname\{Re\}\\operatorname\{tr\}\(Q^\{\\dagger\}Y\)\.The first two terms do not depend onQQ, so the projection problem is equivalent to maximizingRe⁡tr⁡\(Q†​Y\)\\operatorname\{Re\}\\operatorname\{tr\}\(Q^\{\\dagger\}Y\)\. To see why the polar factor attains this maximum, letY=A​Σ​B†Y=A\\Sigma B^\{\\dagger\}be a thin singular value decomposition\. ThenY​\(Y†​Y\)−1/2=A​B†Y\(Y^\{\\dagger\}Y\)^\{\-1/2\}=AB^\{\\dagger\}\. For anyQ∈St​\(p,n\)Q\\in\\mathrm\{St\}\(p,n\),

Re⁡tr⁡\(Q†​Y\)=Re⁡tr⁡\(B†​Q†​A​Σ\)≤∑j=1nσj​\(Y\),\\operatorname\{Re\}\\operatorname\{tr\}\(Q^\{\\dagger\}Y\)=\\operatorname\{Re\}\\operatorname\{tr\}\(B^\{\\dagger\}Q^\{\\dagger\}A\\Sigma\)\\leq\\sum\_\{j=1\}^\{n\}\\sigma\_\{j\}\(Y\),becauseB†​Q†​AB^\{\\dagger\}Q^\{\\dagger\}Ais a contraction and the trace is maximized when its diagonal entries are all equal to one\. This upper bound is achieved by choosingQ=A​B†=R​\(Y\)Q=AB^\{\\dagger\}=R\(Y\), since thenRe⁡tr⁡\(Q†​Y\)=tr⁡\(Σ\)\\operatorname\{Re\}\\operatorname\{tr\}\(Q^\{\\dagger\}Y\)=\\operatorname\{tr\}\(\\Sigma\)\. Thus the polar factorR​\(Y\)=Y​\(Y†​Y\)−1/2R\(Y\)=Y\(Y^\{\\dagger\}Y\)^\{\-1/2\}is the closest Stiefel matrix toYYin Frobenius norm\. Moreover,

R​\(Y\)†​R​\(Y\)=\(Y†​Y\)−1/2​Y†​Y​\(Y†​Y\)−1/2=𝕀n,R\(Y\)^\{\\dagger\}R\(Y\)=\(Y^\{\\dagger\}Y\)^\{\-1/2\}Y^\{\\dagger\}Y\(Y^\{\\dagger\}Y\)^\{\-1/2\}=\\mathbb\{I\}\_\{n\},which is the definition of Stiefel manifold\. ∎

### A\.2Proof of Theorem[1](https://arxiv.org/html/2608.06554#Thmtheorem1)

The following four lemmas, Lemmas[2](https://arxiv.org/html/2608.06554#Thmlemma2)–[5](https://arxiv.org/html/2608.06554#Thmlemma5), will be used to prove Theorem[1](https://arxiv.org/html/2608.06554#Thmtheorem1)\.

###### Lemma 2\.

Suppose Assumption[1](https://arxiv.org/html/2608.06554#Thmassumption1)holds\. ForΓ∈St​\(p,n\)\\Gamma\\in\\mathrm\{St\}\(p,n\), define

SΓ=Γ†​∇F​\(Γ\)\+∇F​\(Γ\)†​Γ2\.S\_\{\\Gamma\}=\\frac\{\\Gamma^\{\\dagger\}\\nabla F\(\\Gamma\)\+\\nabla F\(\\Gamma\)^\{\\dagger\}\\Gamma\}\{2\}\.Then the Riemannian gradient satisfies

grad⁡F​\(Γ\)=∇F​\(Γ\)−Γ​SΓ\.\\operatorname\{grad\}F\(\\Gamma\)=\\nabla F\(\\Gamma\)\-\\Gamma S\_\{\\Gamma\}\.Moreover, for anyΓ,Γ′∈St​\(p,n\)\\Gamma,\\Gamma^\{\\prime\}\\in\\mathrm\{St\}\(p,n\),

F​\(Γ′\)≤\\displaystyle F\(\\Gamma^\{\\prime\}\)\\leqF​\(Γ\)\+⟨grad⁡F​\(Γ\),Γ′−Γ⟩F\+L2​‖Γ′−Γ‖∞2−12​⟨SΓ,\(Γ′−Γ\)†​\(Γ′−Γ\)⟩F\.\\displaystyle F\(\\Gamma\)\+\\left\\langle\\operatorname\{grad\}F\(\\Gamma\),\\Gamma^\{\\prime\}\-\\Gamma\\right\\rangle\_\{F\}\+\\frac\{L\}\{2\}\\\|\\Gamma^\{\\prime\}\-\\Gamma\\\|\_\{\\infty\}^\{2\}\-\\frac\{1\}\{2\}\\left\\langle S\_\{\\Gamma\},\(\\Gamma^\{\\prime\}\-\\Gamma\)^\{\\dagger\}\(\\Gamma^\{\\prime\}\-\\Gamma\)\\right\\rangle\_\{F\}\.Consequently,

F​\(Γ′\)≤F​\(Γ\)\+⟨grad⁡F​\(Γ\),Γ′−Γ⟩F\+L\+‖SΓ‖12​‖Γ′−Γ‖∞2\.F\(\\Gamma^\{\\prime\}\)\\leq F\(\\Gamma\)\+\\left\\langle\\operatorname\{grad\}F\(\\Gamma\),\\Gamma^\{\\prime\}\-\\Gamma\\right\\rangle\_\{F\}\+\\frac\{L\+\\\|S\_\{\\Gamma\}\\\|\_\{1\}\}\{2\}\\\|\\Gamma^\{\\prime\}\-\\Gamma\\\|\_\{\\infty\}^\{2\}\.LettingG1≥supΓ∈St​\(p,n\)‖∇F​\(Γ\)‖1G\_\{1\}\\geq\\sup\_\{\\Gamma\\in\\mathrm\{St\}\(p,n\)\}\\\|\\nabla F\(\\Gamma\)\\\|\_\{1\}, then

F​\(Γ′\)≤F​\(Γ\)\+⟨grad⁡F​\(Γ\),Γ′−Γ⟩F\+L\+G12​‖Γ′−Γ‖∞2\.F\(\\Gamma^\{\\prime\}\)\\leq F\(\\Gamma\)\+\\left\\langle\\operatorname\{grad\}F\(\\Gamma\),\\Gamma^\{\\prime\}\-\\Gamma\\right\\rangle\_\{F\}\+\\frac\{L\+G\_\{1\}\}\{2\}\\\|\\Gamma^\{\\prime\}\-\\Gamma\\\|\_\{\\infty\}^\{2\}\.

###### Proof\.

LetΔ=Γ′−Γ\\Delta=\\Gamma^\{\\prime\}\-\\Gamma\. DefineΓt=Γ\+t​Δ\\Gamma\_\{t\}=\\Gamma\+t\\Delta, for0≤t≤10\\leq t\\leq 1\. By the fundamental theorem of calculus,

F​\(Γ′\)−F​\(Γ\)=∫01⟨∇F​\(Γt\),Δ⟩F​𝑑t\.F\(\\Gamma^\{\\prime\}\)\-F\(\\Gamma\)=\\int\_\{0\}^\{1\}\\left\\langle\\nabla F\(\\Gamma\_\{t\}\),\\Delta\\right\\rangle\_\{F\}\\,dt\.Adding and subtracting∇F​\(Γ\)\\nabla F\(\\Gamma\)inside the integrand gives

F​\(Γ′\)−F​\(Γ\)=⟨∇F​\(Γ\),Δ⟩F\+∫01⟨∇F​\(Γt\)−∇F​\(Γ\),Δ⟩F​𝑑t\.\\displaystyle F\(\\Gamma^\{\\prime\}\)\-F\(\\Gamma\)=\\left\\langle\\nabla F\(\\Gamma\),\\Delta\\right\\rangle\_\{F\}\+\\int\_\{0\}^\{1\}\\left\\langle\\nabla F\(\\Gamma\_\{t\}\)\-\\nabla F\(\\Gamma\),\\Delta\\right\\rangle\_\{F\}\\,dt\.Using spectral–nuclear duality,

⟨∇F​\(Γt\)−∇F​\(Γ\),Δ⟩F≤‖∇F​\(Γt\)−∇F​\(Γ\)‖1​‖Δ‖∞\.\\left\\langle\\nabla F\(\\Gamma\_\{t\}\)\-\\nabla F\(\\Gamma\),\\Delta\\right\\rangle\_\{F\}\\leq\\\|\\nabla F\(\\Gamma\_\{t\}\)\-\\nabla F\(\\Gamma\)\\\|\_\{1\}\\\|\\Delta\\\|\_\{\\infty\}\.By Assumption[1](https://arxiv.org/html/2608.06554#Thmassumption1), applied along the segment joiningΓ\\GammaandΓ′\\Gamma^\{\\prime\},

‖∇F​\(Γt\)−∇F​\(Γ\)‖1≤L​‖Γt−Γ‖∞=L​t​‖Δ‖∞\.\\\|\\nabla F\(\\Gamma\_\{t\}\)\-\\nabla F\(\\Gamma\)\\\|\_\{1\}\\leq L\\\|\\Gamma\_\{t\}\-\\Gamma\\\|\_\{\\infty\}=Lt\\\|\\Delta\\\|\_\{\\infty\}\.Therefore,

∫01⟨∇F​\(Γt\)−∇F​\(Γ\),Δ⟩F​𝑑t≤L​‖Δ‖∞2​∫01t​𝑑t=L2​‖Δ‖∞2\.\\int\_\{0\}^\{1\}\\left\\langle\\nabla F\(\\Gamma\_\{t\}\)\-\\nabla F\(\\Gamma\),\\Delta\\right\\rangle\_\{F\}\\,dt\\leq L\\\|\\Delta\\\|\_\{\\infty\}^\{2\}\\int\_\{0\}^\{1\}t\\,dt=\\frac\{L\}\{2\}\\\|\\Delta\\\|\_\{\\infty\}^\{2\}\.Combining the preceding displays yields

F​\(Γ′\)≤F​\(Γ\)\+⟨∇F​\(Γ\),Δ⟩F\+L2​‖Δ‖∞2\.F\(\\Gamma^\{\\prime\}\)\\leq F\(\\Gamma\)\+\\langle\\nabla F\(\\Gamma\),\\Delta\\rangle\_\{F\}\+\\frac\{L\}\{2\}\\\|\\Delta\\\|\_\{\\infty\}^\{2\}\.The tangent space of the complex Stiefel manifold is

TΓ​St​\(p,n\)=\{ξ∈ℂp×n:Γ†​ξ\+ξ†​Γ=0\},T\_\{\\Gamma\}\\mathrm\{St\}\(p,n\)=\\\{\\xi\\in\\mathbb\{C\}^\{p\\times n\}:\\Gamma^\{\\dagger\}\\xi\+\\xi^\{\\dagger\}\\Gamma=0\\\},which follows by differentiating the constraintΓ†​Γ=𝕀n\\Gamma^\{\\dagger\}\\Gamma=\\mathbb\{I\}\_\{n\}along a differentiable curve onSt​\(p,n\)\\mathrm\{St\}\(p,n\)\. More explicitly, ifγ​\(t\)∈St​\(p,n\)\\gamma\(t\)\\in\\mathrm\{St\}\(p,n\),γ​\(0\)=Γ\\gamma\(0\)=\\Gamma, andγ′​\(0\)=ξ\\gamma^\{\\prime\}\(0\)=\\xi, then differentiatingγ​\(t\)†​γ​\(t\)=𝕀n\\gamma\(t\)^\{\\dagger\}\\gamma\(t\)=\\mathbb\{I\}\_\{n\}att=0t=0givesξ†​Γ\+Γ†​ξ=0\\xi^\{\\dagger\}\\Gamma\+\\Gamma^\{\\dagger\}\\xi=0\.

We next identify the tangent and normal parts of the Euclidean gradient directly\. The matrixSΓS\_\{\\Gamma\}is Hermitian by construction\. Moreover, for any Hermitian matrixSSand anyξ∈TΓ​St​\(p,n\)\\xi\\in T\_\{\\Gamma\}\\mathrm\{St\}\(p,n\),

⟨Γ​S,ξ⟩F=Re⁡tr⁡\(\(Γ​S\)†​ξ\)=Re⁡tr⁡\(S​Γ†​ξ\)\.\\langle\\Gamma S,\\xi\\rangle\_\{F\}=\\operatorname\{Re\}\\operatorname\{tr\}\(\(\\Gamma S\)^\{\\dagger\}\\xi\)=\\operatorname\{Re\}\\operatorname\{tr\}\(S\\Gamma^\{\\dagger\}\\xi\)\.The tangent\-space relation implies thatΓ†​ξ\\Gamma^\{\\dagger\}\\xiis skew\-Hermitian because

\(Γ†​ξ\)†=ξ†​Γ=−Γ†​ξ\.\(\\Gamma^\{\\dagger\}\\xi\)^\{\\dagger\}=\\xi^\{\\dagger\}\\Gamma=\-\\Gamma^\{\\dagger\}\\xi\.For any Hermitian matrixH=H†H=H^\{\\dagger\}and skew\-Hermitian matrixK†=−KK^\{\\dagger\}=\-K, the real Frobenius inner product is zero:

⟨H,K⟩F=Re⁡tr⁡\(H†​K\)=Re⁡tr⁡\(H​K\)\.\\langle H,K\\rangle\_\{F\}=\\operatorname\{Re\}\\operatorname\{tr\}\(H^\{\\dagger\}K\)=\\operatorname\{Re\}\\operatorname\{tr\}\(HK\)\.Since\(tr⁡\(H​K\)\)∗=tr⁡\(\(H​K\)†\)=tr⁡\(K†​H\)=−tr⁡\(K​H\)=−tr⁡\(H​K\)\(\\operatorname\{tr\}\(HK\)\)^\{\\ast\}=\\operatorname\{tr\}\(\(HK\)^\{\\dagger\}\)=\\operatorname\{tr\}\(K^\{\\dagger\}H\)=\-\\operatorname\{tr\}\(KH\)=\-\\operatorname\{tr\}\(HK\), the scalartr⁡\(H​K\)\\operatorname\{tr\}\(HK\)is purely imaginary, and its real part is zero\. Thus⟨Γ​S,ξ⟩F=0\\langle\\Gamma S,\\xi\\rangle\_\{F\}=0for every tangent vectorξ\\xi\. In particular,Γ​SΓ\\Gamma S\_\{\\Gamma\}is normal toTΓ​St​\(p,n\)T\_\{\\Gamma\}\\mathrm\{St\}\(p,n\)\.

It remains to check that the remaining term is tangent\. LetζΓ=∇F​\(Γ\)−Γ​SΓ\\zeta\_\{\\Gamma\}=\\nabla F\(\\Gamma\)\-\\Gamma S\_\{\\Gamma\}\. UsingΓ†​Γ=𝕀n\\Gamma^\{\\dagger\}\\Gamma=\\mathbb\{I\}\_\{n\}andSΓ=SΓ†S\_\{\\Gamma\}=S\_\{\\Gamma\}^\{\\dagger\}, we have

Γ†​ζΓ\+ζΓ†​Γ=Γ†​∇F​\(Γ\)\+∇F​\(Γ\)†​Γ−2​SΓ=0,\\Gamma^\{\\dagger\}\\zeta\_\{\\Gamma\}\+\\zeta\_\{\\Gamma\}^\{\\dagger\}\\Gamma=\\Gamma^\{\\dagger\}\\nabla F\(\\Gamma\)\+\\nabla F\(\\Gamma\)^\{\\dagger\}\\Gamma\-2S\_\{\\Gamma\}=0,where the last equality follows from the definition ofSΓS\_\{\\Gamma\}\. HenceζΓ∈TΓ​St​\(p,n\)\\zeta\_\{\\Gamma\}\\in T\_\{\\Gamma\}\\mathrm\{St\}\(p,n\)\. Therefore

∇F​\(Γ\)=ζΓ\+Γ​SΓ\\nabla F\(\\Gamma\)=\\zeta\_\{\\Gamma\}\+\\Gamma S\_\{\\Gamma\}is an orthogonal decomposition into a tangent component and a normal component\. By the definition of the Riemannian gradient as the tangent projection of the Euclidean gradient,

grad⁡F​\(Γ\)=ΠΓ​\(∇F​\(Γ\)\)=∇F​\(Γ\)−Γ​SΓ\.\\operatorname\{grad\}F\(\\Gamma\)=\\Pi\_\{\\Gamma\}\(\\nabla F\(\\Gamma\)\)=\\nabla F\(\\Gamma\)\-\\Gamma S\_\{\\Gamma\}\.Consequently,

∇F​\(Γ\)=grad⁡F​\(Γ\)\+Γ​SΓ\.\\nabla F\(\\Gamma\)=\\operatorname\{grad\}F\(\\Gamma\)\+\\Gamma S\_\{\\Gamma\}\.Substitution into the Euclidean descent inequality yields

F​\(Γ′\)≤F​\(Γ\)\+⟨grad⁡F​\(Γ\),Δ⟩F\+⟨Γ​SΓ,Δ⟩F\+L2​‖Δ‖∞2\.F\(\\Gamma^\{\\prime\}\)\\leq F\(\\Gamma\)\+\\langle\\operatorname\{grad\}F\(\\Gamma\),\\Delta\\rangle\_\{F\}\+\\langle\\Gamma S\_\{\\Gamma\},\\Delta\\rangle\_\{F\}\+\\frac\{L\}\{2\}\\\|\\Delta\\\|\_\{\\infty\}^\{2\}\.
Now, we are going to express the normal\-gradient contribution as a second\-order term\. Since both endpoints are feasible andΓ′=Γ\+Δ\\Gamma^\{\\prime\}=\\Gamma\+\\Delta,

𝕀n=Γ′⁣†​Γ′=\(Γ\+Δ\)†​\(Γ\+Δ\)=𝕀n\+Γ†​Δ\+Δ†​Γ\+Δ†​Δ,\\mathbb\{I\}\_\{n\}=\\Gamma^\{\\prime\\dagger\}\\Gamma^\{\\prime\}=\(\\Gamma\+\\Delta\)^\{\\dagger\}\(\\Gamma\+\\Delta\)=\\mathbb\{I\}\_\{n\}\+\\Gamma^\{\\dagger\}\\Delta\+\\Delta^\{\\dagger\}\\Gamma\+\\Delta^\{\\dagger\}\\Delta,so

Γ†​Δ\+Δ†​Γ=−Δ†​Δ\.\\Gamma^\{\\dagger\}\\Delta\+\\Delta^\{\\dagger\}\\Gamma=\-\\Delta^\{\\dagger\}\\Delta\.DecomposeΓ†​Δ\\Gamma^\{\\dagger\}\\Deltainto Hermitian and skew\-Hermitian parts:

Γ†​Δ=Γ†​Δ\+Δ†​Γ2\+Γ†​Δ−Δ†​Γ2\.\\Gamma^\{\\dagger\}\\Delta=\\frac\{\\Gamma^\{\\dagger\}\\Delta\+\\Delta^\{\\dagger\}\\Gamma\}\{2\}\+\\frac\{\\Gamma^\{\\dagger\}\\Delta\-\\Delta^\{\\dagger\}\\Gamma\}\{2\}\.The second term is skew\-Hermitian\. SinceSΓS\_\{\\Gamma\}is Hermitian, it is orthogonal to this skew\-Hermitian part under⟨A,B⟩F=Re⁡tr⁡\(A†​B\)\\langle A,B\\rangle\_\{F\}=\\operatorname\{Re\}\\operatorname\{tr\}\(A^\{\\dagger\}B\)\. Therefore, only the Hermitian part ofΓ†​Δ\\Gamma^\{\\dagger\}\\Deltacontributes to the inner product withSΓS\_\{\\Gamma\}\. Using the feasibility identity above,

Γ†​Δ\+Δ†​Γ2=−12​Δ†​Δ\.\\frac\{\\Gamma^\{\\dagger\}\\Delta\+\\Delta^\{\\dagger\}\\Gamma\}\{2\}=\-\\frac\{1\}\{2\}\\Delta^\{\\dagger\}\\Delta\.Also,

⟨Γ​SΓ,Δ⟩F=⟨SΓ,Γ†​Δ⟩F=⟨SΓ,Γ†​Δ\+Δ†​Γ2⟩F=−12​⟨SΓ,Δ†​Δ⟩F\.\\langle\\Gamma S\_\{\\Gamma\},\\Delta\\rangle\_\{F\}=\\langle S\_\{\\Gamma\},\\Gamma^\{\\dagger\}\\Delta\\rangle\_\{F\}=\\left\\langle S\_\{\\Gamma\},\\frac\{\\Gamma^\{\\dagger\}\\Delta\+\\Delta^\{\\dagger\}\\Gamma\}\{2\}\\right\\rangle\_\{F\}=\-\\frac\{1\}\{2\}\\langle S\_\{\\Gamma\},\\Delta^\{\\dagger\}\\Delta\\rangle\_\{F\}\.Substituting this identity into the descent inequality proves the first bound\. The correction is quadratic because it depends onΔ†​Δ\\Delta^\{\\dagger\}\\Delta, even thoughΔ=Γ′−Γ\\Delta=\\Gamma^\{\\prime\}\-\\Gammais generally not a tangent vector atΓ\\Gamma\.

For the second bound, use spectral–nuclear duality:

−12​⟨SΓ,Δ†​Δ⟩F≤12​\|⟨SΓ,Δ†​Δ⟩F\|≤12​‖SΓ‖1​‖Δ†​Δ‖∞\.\-\\frac\{1\}\{2\}\\langle S\_\{\\Gamma\},\\Delta^\{\\dagger\}\\Delta\\rangle\_\{F\}\\leq\\frac\{1\}\{2\}\\left\|\\langle S\_\{\\Gamma\},\\Delta^\{\\dagger\}\\Delta\\rangle\_\{F\}\\right\|\\leq\\frac\{1\}\{2\}\\\|S\_\{\\Gamma\}\\\|\_\{1\}\\\|\\Delta^\{\\dagger\}\\Delta\\\|\_\{\\infty\}\.Furthermore, by submultiplicativity,

‖Δ†​Δ‖∞≤‖Δ†‖∞​‖Δ‖∞=‖Δ‖∞2\.\\\|\\Delta^\{\\dagger\}\\Delta\\\|\_\{\\infty\}\\leq\\\|\\Delta^\{\\dagger\}\\\|\_\{\\infty\}\\\|\\Delta\\\|\_\{\\infty\}=\\\|\\Delta\\\|\_\{\\infty\}^\{2\}\.Combining these two estimates gives the second bound\.

Finally, we control‖SΓ‖1\\\|S\_\{\\Gamma\}\\\|\_\{1\}\. SinceSΓ=\(Γ†​∇F​\(Γ\)\+∇F​\(Γ\)†​Γ\)/2S\_\{\\Gamma\}=\(\\Gamma^\{\\dagger\}\\nabla F\(\\Gamma\)\+\\nabla F\(\\Gamma\)^\{\\dagger\}\\Gamma\)/2, the triangle inequality gives

‖SΓ‖1≤12​‖Γ†​∇F​\(Γ\)‖1\+12​‖∇F​\(Γ\)†​Γ‖1\.\\\|S\_\{\\Gamma\}\\\|\_\{1\}\\leq\\frac\{1\}\{2\}\\\|\\Gamma^\{\\dagger\}\\nabla F\(\\Gamma\)\\\|\_\{1\}\+\\frac\{1\}\{2\}\\\|\\nabla F\(\\Gamma\)^\{\\dagger\}\\Gamma\\\|\_\{1\}\.The two terms are equal, and becauseΓ†​Γ=𝕀n\\Gamma^\{\\dagger\}\\Gamma=\\mathbb\{I\}\_\{n\}, all singular values ofΓ\\Gammaare equal to one; hence‖Γ†‖∞=1\\\|\\Gamma^\{\\dagger\}\\\|\_\{\\infty\}=1\. Therefore,

‖SΓ‖1≤‖Γ†​∇F​\(Γ\)‖1≤‖∇F​\(Γ\)‖1≤G1\.\\\|S\_\{\\Gamma\}\\\|\_\{1\}\\leq\\\|\\Gamma^\{\\dagger\}\\nabla F\(\\Gamma\)\\\|\_\{1\}\\leq\\\|\\nabla F\(\\Gamma\)\\\|\_\{1\}\\leq G\_\{1\}\.Substituting this estimate into the second bound gives the final statement\. ∎

The following Lemma[3](https://arxiv.org/html/2608.06554#Thmlemma3)shows that the Euclidean Lipschitz smoothness assumption implies a corresponding Lipschitz bound for the Riemannian gradient after projection onto the Stiefel tangent space\.

###### Lemma 3\.

Suppose Assumption[1](https://arxiv.org/html/2608.06554#Thmassumption1)holds, and letG1≥supΓ∈St​\(p,n\)‖∇F​\(Γ\)‖1G\_\{1\}\\geq\\sup\_\{\\Gamma\\in\\mathrm\{St\}\(p,n\)\}\\\|\\nabla F\(\\Gamma\)\\\|\_\{1\}as in Lemma[2](https://arxiv.org/html/2608.06554#Thmlemma2)\. Then, for allΓ,Γ′∈St​\(p,n\)\\Gamma,\\Gamma^\{\\prime\}\\in\\mathrm\{St\}\(p,n\),

‖grad⁡F​\(Γ\)−grad⁡F​\(Γ′\)‖1≤LR​‖Γ−Γ′‖∞,LR:=2​\(L\+G1\)\.\\\|\\operatorname\{grad\}F\(\\Gamma\)\-\\operatorname\{grad\}F\(\\Gamma^\{\\prime\}\)\\\|\_\{1\}\\leq L\_\{R\}\\\|\\Gamma\-\\Gamma^\{\\prime\}\\\|\_\{\\infty\},\\qquad L\_\{R\}:=2\(L\+G\_\{1\}\)\.

###### Proof\.

LetA=∇F​\(Γ\)A=\\nabla F\(\\Gamma\),A′=∇F​\(Γ′\)A^\{\\prime\}=\\nabla F\(\\Gamma^\{\\prime\}\), andd=‖Γ−Γ′‖∞d=\\\|\\Gamma\-\\Gamma^\{\\prime\}\\\|\_\{\\infty\}\. By Assumption[1](https://arxiv.org/html/2608.06554#Thmassumption1),

‖A−A′‖1≤L​d\.\\\|A\-A^\{\\prime\}\\\|\_\{1\}\\leq Ld\.Usinggrad⁡F​\(Γ\)=∇F​\(Γ\)−Γ​SΓ\\operatorname\{grad\}F\(\\Gamma\)=\\nabla F\(\\Gamma\)\-\\Gamma S\_\{\\Gamma\}, we have

‖grad⁡F​\(Γ\)−grad⁡F​\(Γ′\)‖1\\displaystyle\\\|\\operatorname\{grad\}F\(\\Gamma\)\-\\operatorname\{grad\}F\(\\Gamma^\{\\prime\}\)\\\|\_\{1\}≤‖A−A′‖1\+‖Γ​SΓ−Γ′​SΓ′‖1\\displaystyle\\leq\\\|A\-A^\{\\prime\}\\\|\_\{1\}\+\\\|\\Gamma S\_\{\\Gamma\}\-\\Gamma^\{\\prime\}S\_\{\\Gamma^\{\\prime\}\}\\\|\_\{1\}≤L​d\+‖\(Γ−Γ′\)​SΓ‖1\+‖Γ′​\(SΓ−SΓ′\)‖1\.\\displaystyle\\leq Ld\+\\\|\(\\Gamma\-\\Gamma^\{\\prime\}\)S\_\{\\Gamma\}\\\|\_\{1\}\+\\\|\\Gamma^\{\\prime\}\(S\_\{\\Gamma\}\-S\_\{\\Gamma^\{\\prime\}\}\)\\\|\_\{1\}\.SinceΓ,Γ′∈St​\(p,n\)\\Gamma,\\Gamma^\{\\prime\}\\in\\mathrm\{St\}\(p,n\),‖Γ′‖∞=1\\\|\\Gamma^\{\\prime\}\\\|\_\{\\infty\}=1\. Also‖SΓ‖1≤‖∇F​\(Γ\)‖1≤G1\\\|S\_\{\\Gamma\}\\\|\_\{1\}\\leq\\\|\\nabla F\(\\Gamma\)\\\|\_\{1\}\\leq G\_\{1\}, as shown in Lemma[2](https://arxiv.org/html/2608.06554#Thmlemma2)\. Therefore, by submultiplicativity of the nuclear norm with respect to the operator norm,

‖\(Γ−Γ′\)​SΓ‖1≤‖Γ−Γ′‖∞​‖SΓ‖1≤G1​d,\\\|\(\\Gamma\-\\Gamma^\{\\prime\}\)S\_\{\\Gamma\}\\\|\_\{1\}\\leq\\\|\\Gamma\-\\Gamma^\{\\prime\}\\\|\_\{\\infty\}\\\|S\_\{\\Gamma\}\\\|\_\{1\}\\leq G\_\{1\}d,and

‖Γ′​\(SΓ−SΓ′\)‖1≤‖Γ′‖∞​‖SΓ−SΓ′‖1=‖SΓ−SΓ′‖1\.\\\|\\Gamma^\{\\prime\}\(S\_\{\\Gamma\}\-S\_\{\\Gamma^\{\\prime\}\}\)\\\|\_\{1\}\\leq\\\|\\Gamma^\{\\prime\}\\\|\_\{\\infty\}\\\|S\_\{\\Gamma\}\-S\_\{\\Gamma^\{\\prime\}\}\\\|\_\{1\}=\\\|S\_\{\\Gamma\}\-S\_\{\\Gamma^\{\\prime\}\}\\\|\_\{1\}\.Combining these two estimates with‖A−A′‖1≤L​d\\\|A\-A^\{\\prime\}\\\|\_\{1\}\\leq Ldgives

‖grad⁡F​\(Γ\)−grad⁡F​\(Γ′\)‖1≤\(L\+G1\)​d\+‖SΓ−SΓ′‖1\.\\\|\\operatorname\{grad\}F\(\\Gamma\)\-\\operatorname\{grad\}F\(\\Gamma^\{\\prime\}\)\\\|\_\{1\}\\leq\(L\+G\_\{1\}\)d\+\\\|S\_\{\\Gamma\}\-S\_\{\\Gamma^\{\\prime\}\}\\\|\_\{1\}\.It remains to bound the last term\. From the definition ofSΓS\_\{\\Gamma\}, with

B=Γ†​A−Γ′⁣†​A′,B=\\Gamma^\{\\dagger\}A\-\\Gamma^\{\\prime\\dagger\}A^\{\\prime\},we have

SΓ−SΓ′=12​\(B\+B†\),S\_\{\\Gamma\}\-S\_\{\\Gamma^\{\\prime\}\}=\\frac\{1\}\{2\}\\left\(B\+B^\{\\dagger\}\\right\),because

B†=A†​Γ−A′⁣†​Γ′\.B^\{\\dagger\}=A^\{\\dagger\}\\Gamma\-A^\{\\prime\\dagger\}\\Gamma^\{\\prime\}\.Since the nuclear norm is invariant under adjoints,‖B†‖1=‖B‖1\\\|B^\{\\dagger\}\\\|\_\{1\}=\\\|B\\\|\_\{1\}\. Thus

‖SΓ−SΓ′‖1≤12​‖B‖1\+12​‖B†‖1=‖B‖1\.\\\|S\_\{\\Gamma\}\-S\_\{\\Gamma^\{\\prime\}\}\\\|\_\{1\}\\leq\\frac\{1\}\{2\}\\\|B\\\|\_\{1\}\+\\frac\{1\}\{2\}\\\|B^\{\\dagger\}\\\|\_\{1\}=\\\|B\\\|\_\{1\}\.That is,

‖SΓ−SΓ′‖1\\displaystyle\\\|S\_\{\\Gamma\}\-S\_\{\\Gamma^\{\\prime\}\}\\\|\_\{1\}≤‖Γ†​A−Γ′⁣†​A′‖1\\displaystyle\\leq\\\|\\Gamma^\{\\dagger\}A\-\\Gamma^\{\\prime\\dagger\}A^\{\\prime\}\\\|\_\{1\}≤‖Γ†​\(A−A′\)‖1\+‖\(Γ−Γ′\)†​A′‖1\\displaystyle\\leq\\\|\\Gamma^\{\\dagger\}\(A\-A^\{\\prime\}\)\\\|\_\{1\}\+\\\|\(\\Gamma\-\\Gamma^\{\\prime\}\)^\{\\dagger\}A^\{\\prime\}\\\|\_\{1\}≤L​d\+G1​d\.\\displaystyle\\leq Ld\+G\_\{1\}d\.Combining the two estimates gives

‖grad⁡F​\(Γ\)−grad⁡F​\(Γ′\)‖1≤2​\(L\+G1\)​d=LR​‖Γ−Γ′‖∞\.\\\|\\operatorname\{grad\}F\(\\Gamma\)\-\\operatorname\{grad\}F\(\\Gamma^\{\\prime\}\)\\\|\_\{1\}\\leq 2\(L\+G\_\{1\}\)d=L\_\{R\}\\\|\\Gamma\-\\Gamma^\{\\prime\}\\\|\_\{\\infty\}\.∎

The following Lemma[4](https://arxiv.org/html/2608.06554#Thmlemma4)gives the one\-step descent estimate for the actual NS\-RIS update, including both the finite Newton–Schulz direction error and the second Newton–Schulz feasibility correction\.

###### Lemma 4\.

Consider the actual NS\-RIS update

Γk\+12=Γk−η​M~k\+1,Γk\+1=NSOrth⁡\(Γk\+12,TNS\),\\Gamma\_\{k\+\\frac\{1\}\{2\}\}=\\Gamma\_\{k\}\-\\eta\\widetilde\{M\}\_\{k\+1\},\\qquad\\Gamma\_\{k\+1\}=\\operatorname\{NSOrth\}\(\\Gamma\_\{k\+\\frac\{1\}\{2\}\},T\_\{\\mathrm\{NS\}\}\),and suppose thatΓk\+1∈St​\(p,n\)\\Gamma\_\{k\+1\}\\in\\mathrm\{St\}\(p,n\)\. Define the second Newton–Schulz correction

Ek\+1=Γk\+1−Γk\+12,ρk\+1=‖Ek\+1‖∞\.E\_\{k\+1\}=\\Gamma\_\{k\+1\}\-\\Gamma\_\{k\+\\frac\{1\}\{2\}\},\\qquad\\rho\_\{k\+1\}=\\\|E\_\{k\+1\}\\\|\_\{\\infty\}\.Then, under Assumptions[1](https://arxiv.org/html/2608.06554#Thmassumption1)and[3](https://arxiv.org/html/2608.06554#Thmassumption3),

F​\(Γk\+1\)≤\\displaystyle F\(\\Gamma\_\{k\+1\}\)\\leqF​\(Γk\)−η​\(1−εNS\)​‖Mk\+1‖1\\displaystyle F\(\\Gamma\_\{k\}\)\-\\eta\(1\-\\varepsilon\_\{\\mathrm\{NS\}\}\)\\\|M\_\{k\+1\}\\\|\_\{1\}\+η​\(1\+εNS\)​‖grad⁡F​\(Γk\)−Mk\+1‖1\+ρk\+1​‖grad⁡F​\(Γk\)‖1\\displaystyle\+\\eta\(1\+\\varepsilon\_\{\\mathrm\{NS\}\}\)\\\|\\operatorname\{grad\}F\(\\Gamma\_\{k\}\)\-M\_\{k\+1\}\\\|\_\{1\}\+\\rho\_\{k\+1\}\\\|\\operatorname\{grad\}F\(\\Gamma\_\{k\}\)\\\|\_\{1\}\+LR2​\(η​‖M~k\+1‖∞\+ρk\+1\)2\.\\displaystyle\+\\frac\{L\_\{R\}\}\{2\}\\left\(\\eta\\\|\\widetilde\{M\}\_\{k\+1\}\\\|\_\{\\infty\}\+\\rho\_\{k\+1\}\\right\)^\{2\}\.

###### Proof\.

Let

Hk=grad⁡F​\(Γk\)andΔk=Γk\+1−Γk\.H\_\{k\}=\\operatorname\{grad\}F\(\\Gamma\_\{k\}\)\\quad\\text\{and\}\\quad\\Delta\_\{k\}=\\Gamma\_\{k\+1\}\-\\Gamma\_\{k\}\.The actual update gives the exact decomposition

Δk=−η​M~k\+1\+Ek\+1\.\\Delta\_\{k\}=\-\\eta\\widetilde\{M\}\_\{k\+1\}\+E\_\{k\+1\}\.By Lemma[2](https://arxiv.org/html/2608.06554#Thmlemma2), with the weaker quadratic constantLR≥L\+G1L\_\{R\}\\geq L\+G\_\{1\},

F​\(Γk\+1\)≤F​\(Γk\)−η​⟨Hk,M~k\+1⟩F\+⟨Hk,Ek\+1⟩F\+LR2​‖Δk‖∞2\.F\(\\Gamma\_\{k\+1\}\)\\leq F\(\\Gamma\_\{k\}\)\-\\eta\\langle H\_\{k\},\\widetilde\{M\}\_\{k\+1\}\\rangle\_\{F\}\+\\langle H\_\{k\},E\_\{k\+1\}\\rangle\_\{F\}\+\\frac\{L\_\{R\}\}\{2\}\\\|\\Delta\_\{k\}\\\|\_\{\\infty\}^\{2\}\.Spectral–nuclear duality and the definition ofρk\+1\\rho\_\{k\+1\}give

⟨Hk,Ek\+1⟩F≤‖Hk‖1​‖Ek\+1‖∞=ρk\+1​‖Hk‖1,\\langle H\_\{k\},E\_\{k\+1\}\\rangle\_\{F\}\\leq\\\|H\_\{k\}\\\|\_\{1\}\\\|E\_\{k\+1\}\\\|\_\{\\infty\}=\\rho\_\{k\+1\}\\\|H\_\{k\}\\\|\_\{1\},and the triangle inequality gives

‖Δk‖∞≤η​‖M~k\+1‖∞\+ρk\+1\.\\\|\\Delta\_\{k\}\\\|\_\{\\infty\}\\leq\\eta\\\|\\widetilde\{M\}\_\{k\+1\}\\\|\_\{\\infty\}\+\\rho\_\{k\+1\}\.The inner product is decomposed into a momentum term and a tracking error:

⟨Hk,M~k\+1⟩F\\displaystyle\\langle H\_\{k\},\\widetilde\{M\}\_\{k\+1\}\\rangle\_\{F\}=⟨Mk\+1,M~k\+1⟩F\+⟨Hk−Mk\+1,M~k\+1⟩F\.\\displaystyle=\\langle M\_\{k\+1\},\\widetilde\{M\}\_\{k\+1\}\\rangle\_\{F\}\+\\langle H\_\{k\}\-M\_\{k\+1\},\\widetilde\{M\}\_\{k\+1\}\\rangle\_\{F\}\.LetPk\+1=Polar⁡\(Mk\+1\)P\_\{k\+1\}=\\operatorname\{Polar\}\(M\_\{k\+1\}\)\. Since⟨Mk\+1,Pk\+1⟩F=‖Mk\+1‖1\\langle M\_\{k\+1\},P\_\{k\+1\}\\rangle\_\{F\}=\\\|M\_\{k\+1\}\\\|\_\{1\},

⟨Mk\+1,M~k\+1⟩F=‖Mk\+1‖1\+⟨Mk\+1,M~k\+1−Pk\+1⟩F≥\(1−εNS\)​‖Mk\+1‖1,\\langle M\_\{k\+1\},\\widetilde\{M\}\_\{k\+1\}\\rangle\_\{F\}=\\\|M\_\{k\+1\}\\\|\_\{1\}\+\\langle M\_\{k\+1\},\\widetilde\{M\}\_\{k\+1\}\-P\_\{k\+1\}\\rangle\_\{F\}\\geq\(1\-\\varepsilon\_\{\\mathrm\{NS\}\}\)\\\|M\_\{k\+1\}\\\|\_\{1\},where the last inequality uses spectral–nuclear duality and Assumption[3](https://arxiv.org/html/2608.06554#Thmassumption3)\. The same assumption also implies

‖M~k\+1‖∞≤‖Pk\+1‖∞\+‖M~k\+1−Pk\+1‖∞≤1\+εNS\.\\\|\\widetilde\{M\}\_\{k\+1\}\\\|\_\{\\infty\}\\leq\\\|P\_\{k\+1\}\\\|\_\{\\infty\}\+\\\|\\widetilde\{M\}\_\{k\+1\}\-P\_\{k\+1\}\\\|\_\{\\infty\}\\leq 1\+\\varepsilon\_\{\\mathrm\{NS\}\}\.Thus the remaining term is bounded by

\|⟨Hk−Mk\+1,M~k\+1⟩F\|≤\(1\+εNS\)​‖Hk−Mk\+1‖1\.\\left\|\\langle H\_\{k\}\-M\_\{k\+1\},\\widetilde\{M\}\_\{k\+1\}\\rangle\_\{F\}\\right\|\\leq\(1\+\\varepsilon\_\{\\mathrm\{NS\}\}\)\\\|H\_\{k\}\-M\_\{k\+1\}\\\|\_\{1\}\.Combining the inequalities proves the claim\. ∎

The following Lemma[5](https://arxiv.org/html/2608.06554#Thmlemma5)controls how closely the momentum variable tracks the Riemannian gradient along the stochastic NS\-RIS iterates\.

###### Lemma 5\.

For a constant stepsizeη\\etaand a fixed momentum parameter0≤β<10\\leq\\beta<1, let

Hk=grad⁡F​\(Γk\),ek\+1=Mk\+1−Hk\.H\_\{k\}=\\operatorname\{grad\}F\(\\Gamma\_\{k\}\),\\qquad e\_\{k\+1\}=M\_\{k\+1\}\-H\_\{k\}\.Fork≥1k\\geq 1, define the second Newton–Schulz correction

Ek=Γk−\(Γk−1−η​M~k\),ρk=‖Ek‖∞\.E\_\{k\}=\\Gamma\_\{k\}\-\\bigl\(\\Gamma\_\{k\-1\}\-\\eta\\widetilde\{M\}\_\{k\}\\bigr\),\\qquad\\rho\_\{k\}=\\\|E\_\{k\}\\\|\_\{\\infty\}\.Under Assumptions[1](https://arxiv.org/html/2608.06554#Thmassumption1)and[2](https://arxiv.org/html/2608.06554#Thmassumption2), the exponential moving averageMk\+1=β​Mk\+\(1−β\)​gkM\_\{k\+1\}=\\beta M\_\{k\}\+\(1\-\\beta\)g\_\{k\}satisfies

1K​∑k=0K−1𝔼​‖Hk−Mk\+1‖1≤\\displaystyle\\frac\{1\}\{K\}\\sum\_\{k=0\}^\{K\-1\}\\mathbb\{E\}\\\|H\_\{k\}\-M\_\{k\+1\}\\\|\_\{1\}\\leq𝔼​‖e1‖1\(1−β\)​K\+n​σb\+β​LR1−β​1K​∑k=1K−1𝔼​\[η​‖M~k‖∞\+ρk\]\.\\displaystyle\\frac\{\\mathbb\{E\}\\\|e\_\{1\}\\\|\_\{1\}\}\{\(1\-\\beta\)K\}\+\\frac\{\\sqrt\{n\}\\,\\sigma\}\{\\sqrt\{b\}\}\+\\frac\{\\beta L\_\{R\}\}\{1\-\\beta\}\\frac\{1\}\{K\}\\sum\_\{k=1\}^\{K\-1\}\\mathbb\{E\}\\left\[\\eta\\\|\\widetilde\{M\}\_\{k\}\\\|\_\{\\infty\}\+\\rho\_\{k\}\\right\]\.

###### Proof\.

Write

gk=Hk\+ζk,where𝔼​\[ζk∣ℱk\]=0\.g\_\{k\}=H\_\{k\}\+\\zeta\_\{k\},\\qquad\\text\{where\}\\quad\\mathbb\{E\}\[\\zeta\_\{k\}\\mid\\mathcal\{F\}\_\{k\}\]=0\.By Assumption[2](https://arxiv.org/html/2608.06554#Thmassumption2),

𝔼​\[‖ζk‖F2∣ℱk\]≤σ2b\.\\mathbb\{E\}\\\!\\left\[\\\|\\zeta\_\{k\}\\\|\_\{F\}^\{2\}\\mid\\mathcal\{F\}\_\{k\}\\right\]\\leq\\frac\{\\sigma^\{2\}\}\{b\}\.Sincerank⁡\(ζk\)≤n\\operatorname\{rank\}\(\\zeta\_\{k\}\)\\leq n,

𝔼​‖ζk‖1≤n​𝔼​‖ζk‖F≤n​σb\.\\mathbb\{E\}\\\|\\zeta\_\{k\}\\\|\_\{1\}\\leq\\sqrt\{n\}\\,\\mathbb\{E\}\\\|\\zeta\_\{k\}\\\|\_\{F\}\\leq\\frac\{\\sqrt\{n\}\\,\\sigma\}\{\\sqrt\{b\}\}\.Fork≥1k\\geq 1, the momentum recursion gives

ek\+1\\displaystyle e\_\{k\+1\}=Mk\+1−Hk\\displaystyle=M\_\{k\+1\}\-H\_\{k\}=β​Mk\+\(1−β\)​gk−Hk\\displaystyle=\\beta M\_\{k\}\+\(1\-\\beta\)g\_\{k\}\-H\_\{k\}=β​\(Mk−Hk−1\)\+β​\(Hk−1−Hk\)\+\(1−β\)​ζk\\displaystyle=\\beta\(M\_\{k\}\-H\_\{k\-1\}\)\+\\beta\(H\_\{k\-1\}\-H\_\{k\}\)\+\(1\-\\beta\)\\zeta\_\{k\}=β​ek\+β​\(Hk−1−Hk\)\+\(1−β\)​ζk\.\\displaystyle=\\beta e\_\{k\}\+\\beta\(H\_\{k\-1\}\-H\_\{k\}\)\+\(1\-\\beta\)\\zeta\_\{k\}\.Hence,

𝔼​‖ek\+1‖1≤β​𝔼​‖ek‖1\+β​𝔼​‖Hk−Hk−1‖1\+\(1−β\)​n​σb\.\\mathbb\{E\}\\\|e\_\{k\+1\}\\\|\_\{1\}\\leq\\beta\\mathbb\{E\}\\\|e\_\{k\}\\\|\_\{1\}\+\\beta\\mathbb\{E\}\\\|H\_\{k\}\-H\_\{k\-1\}\\\|\_\{1\}\+\(1\-\\beta\)\\frac\{\\sqrt\{n\}\\,\\sigma\}\{\\sqrt\{b\}\}\.By Lemma[3](https://arxiv.org/html/2608.06554#Thmlemma3)and the actual update,

Γk−Γk−1=−η​M~k\+Ek,\\Gamma\_\{k\}\-\\Gamma\_\{k\-1\}=\-\\eta\\widetilde\{M\}\_\{k\}\+E\_\{k\},we have

‖Hk−Hk−1‖1≤LR​‖Γk−Γk−1‖∞≤LR​\(η​‖M~k‖∞\+ρk\)\.\\\|H\_\{k\}\-H\_\{k\-1\}\\\|\_\{1\}\\leq L\_\{R\}\\\|\\Gamma\_\{k\}\-\\Gamma\_\{k\-1\}\\\|\_\{\\infty\}\\leq L\_\{R\}\\left\(\\eta\\\|\\widetilde\{M\}\_\{k\}\\\|\_\{\\infty\}\+\\rho\_\{k\}\\right\)\.Letak=𝔼​‖ek‖1a\_\{k\}=\\mathbb\{E\}\\\|e\_\{k\}\\\|\_\{1\}ands=n​σ/bs=\\sqrt\{n\}\\sigma/\\sqrt\{b\}\. Summing the preceding recursion fromk=1k=1toK−1K\-1gives

∑k=1K−1ak\+1≤β​∑k=1K−1ak\+β​LR​∑k=1K−1𝔼​\[η​‖M~k‖∞\+ρk\]\+\(K−1\)​\(1−β\)​s\.\\sum\_\{k=1\}^\{K\-1\}a\_\{k\+1\}\\leq\\beta\\sum\_\{k=1\}^\{K\-1\}a\_\{k\}\+\\beta L\_\{R\}\\sum\_\{k=1\}^\{K\-1\}\\mathbb\{E\}\\left\[\\eta\\\|\\widetilde\{M\}\_\{k\}\\\|\_\{\\infty\}\+\\rho\_\{k\}\\right\]\+\(K\-1\)\(1\-\\beta\)s\.Since∑k=1K−1ak\+1=∑k=2Kak\\sum\_\{k=1\}^\{K\-1\}a\_\{k\+1\}=\\sum\_\{k=2\}^\{K\}a\_\{k\}and∑k=1K−1ak=a1\+∑k=2Kak−aK\\sum\_\{k=1\}^\{K\-1\}a\_\{k\}=a\_\{1\}\+\\sum\_\{k=2\}^\{K\}a\_\{k\}\-a\_\{K\}, this implies

\(1−β\)​∑k=2Kak≤β​a1\+β​LR​∑k=1K−1𝔼​\[η​‖M~k‖∞\+ρk\]\+\(K−1\)​\(1−β\)​s,\(1\-\\beta\)\\sum\_\{k=2\}^\{K\}a\_\{k\}\\leq\\beta a\_\{1\}\+\\beta L\_\{R\}\\sum\_\{k=1\}^\{K\-1\}\\mathbb\{E\}\\left\[\\eta\\\|\\widetilde\{M\}\_\{k\}\\\|\_\{\\infty\}\+\\rho\_\{k\}\\right\]\+\(K\-1\)\(1\-\\beta\)s,where the nonpositive term−β​aK\-\\beta a\_\{K\}has been dropped\. Addinga1a\_\{1\}to both sides after division by1−β1\-\\betagives

∑k=1Kak≤a11−β\+\(K−1\)​s\+β​LR1−β​∑k=1K−1𝔼​\[η​‖M~k‖∞\+ρk\]\.\\sum\_\{k=1\}^\{K\}a\_\{k\}\\leq\\frac\{a\_\{1\}\}\{1\-\\beta\}\+\(K\-1\)s\+\\frac\{\\beta L\_\{R\}\}\{1\-\\beta\}\\sum\_\{k=1\}^\{K\-1\}\\mathbb\{E\}\\left\[\\eta\\\|\\widetilde\{M\}\_\{k\}\\\|\_\{\\infty\}\+\\rho\_\{k\}\\right\]\.After the final division byKK, the noise contribution is\(\(K−1\)/K\)​s≤s\(\(K\-1\)/K\)s\\leq s\. Dividing byKKand recalling that‖Hk−Mk\+1‖1=‖ek\+1‖1\\\|H\_\{k\}\-M\_\{k\+1\}\\\|\_\{1\}=\\\|e\_\{k\+1\}\\\|\_\{1\}proves the claim\. ∎

###### Proof of Theorem[1](https://arxiv.org/html/2608.06554#Thmtheorem1)\.

For compactness, write

Hk=grad⁡F​\(Γk\)andTk=‖Hk−Mk\+1‖1\.H\_\{k\}=\\operatorname\{grad\}F\(\\Gamma\_\{k\}\)\\quad\\text\{and\}\\quad T\_\{k\}=\\\|H\_\{k\}\-M\_\{k\+1\}\\\|\_\{1\}\.Lemma[4](https://arxiv.org/html/2608.06554#Thmlemma4)gives the one\-step inequality

F​\(Γk\+1\)≤\\displaystyle F\(\\Gamma\_\{k\+1\}\)\\leqF​\(Γk\)−η​\(1−εNS\)​‖Mk\+1‖1\+η​\(1\+εNS\)​Tk\+ρk\+1​‖Hk‖1\+LR2​\(η​‖M~k\+1‖∞\+ρk\+1\)2\.\\displaystyle F\(\\Gamma\_\{k\}\)\-\\eta\(1\-\\varepsilon\_\{\\mathrm\{NS\}\}\)\\\|M\_\{k\+1\}\\\|\_\{1\}\+\\eta\(1\+\\varepsilon\_\{\\mathrm\{NS\}\}\)T\_\{k\}\+\\rho\_\{k\+1\}\\\|H\_\{k\}\\\|\_\{1\}\+\\frac\{L\_\{R\}\}\{2\}\\left\(\\eta\\\|\\widetilde\{M\}\_\{k\+1\}\\\|\_\{\\infty\}\+\\rho\_\{k\+1\}\\right\)^\{2\}\.By Assumption[3](https://arxiv.org/html/2608.06554#Thmassumption3),‖M~k\+1‖∞≤\(1\+εNS\)\\\|\\widetilde\{M\}\_\{k\+1\}\\\|\_\{\\infty\}\\leq\(1\+\\varepsilon\_\{\\mathrm\{NS\}\}\)\. Also,‖Hk‖1≤‖∇F​\(Γk\)‖1\+‖SΓk‖1≤2​G1\\\|H\_\{k\}\\\|\_\{1\}\\leq\\\|\\nabla F\(\\Gamma\_\{k\}\)\\\|\_\{1\}\+\\\|S\_\{\\Gamma\_\{k\}\}\\\|\_\{1\}\\leq 2G\_\{1\}\. Using\(x\+y\)2≤2​x2\+2​y2\(x\+y\)^\{2\}\\leq 2x^\{2\}\+2y^\{2\}gives

η​\(1−εNS\)​‖Mk\+1‖1≤F​\(Γk\)−F​\(Γk\+1\)\+η​\(1\+εNS\)​Tk\+2​G1​ρk\+1\+LR​η2​\(1\+εNS\)2\+LR​ρk\+12\.\\eta\(1\-\\varepsilon\_\{\\mathrm\{NS\}\}\)\\\|M\_\{k\+1\}\\\|\_\{1\}\\leq F\(\\Gamma\_\{k\}\)\-F\(\\Gamma\_\{k\+1\}\)\+\\eta\(1\+\\varepsilon\_\{\\mathrm\{NS\}\}\)T\_\{k\}\+2G\_\{1\}\\rho\_\{k\+1\}\+L\_\{R\}\\eta^\{2\}\(1\+\\varepsilon\_\{\\mathrm\{NS\}\}\)^\{2\}\+L\_\{R\}\\rho\_\{k\+1\}^\{2\}\.Taking expectations, summing fromk=0k=0toK−1K\-1, and usingF​\(ΓK\)≥F∗F\(\\Gamma\_\{K\}\)\\geq F^\{\*\}gives

1K​∑k=0K−1𝔼​‖Mk\+1‖1≤Dη​\(1−εNS\)​K\+1\+εNS1−εNS​1K​∑k=0K−1𝔼​Tk\+2​G1​ρ¯K\+LR​q¯Kη​\(1−εNS\)\+LR​\(1\+εNS\)2​η1−εNS\.\\frac\{1\}\{K\}\\sum\_\{k=0\}^\{K\-1\}\\mathbb\{E\}\\\|M\_\{k\+1\}\\\|\_\{1\}\\leq\\frac\{D\}\{\\eta\(1\-\\varepsilon\_\{\\mathrm\{NS\}\}\)K\}\+\\frac\{1\+\\varepsilon\_\{\\mathrm\{NS\}\}\}\{1\-\\varepsilon\_\{\\mathrm\{NS\}\}\}\\frac\{1\}\{K\}\\sum\_\{k=0\}^\{K\-1\}\\mathbb\{E\}T\_\{k\}\+\\frac\{2G\_\{1\}\\bar\{\\rho\}\_\{K\}\+L\_\{R\}\\bar\{q\}\_\{K\}\}\{\\eta\(1\-\\varepsilon\_\{\\mathrm\{NS\}\}\)\}\+\\frac\{L\_\{R\}\(1\+\\varepsilon\_\{\\mathrm\{NS\}\}\)^\{2\}\\eta\}\{1\-\\varepsilon\_\{\\mathrm\{NS\}\}\}\.The desired stationarity measure is related to the momentum by the triangle inequality:

‖Hk‖1≤‖Mk\+1‖1\+Tk\.\\\|H\_\{k\}\\\|\_\{1\}\\leq\\\|M\_\{k\+1\}\\\|\_\{1\}\+T\_\{k\}\.Combining the last two inequalities gives

1K​∑k=0K−1𝔼​‖Hk‖1≤\\displaystyle\\frac\{1\}\{K\}\\sum\_\{k=0\}^\{K\-1\}\\mathbb\{E\}\\\|H\_\{k\}\\\|\_\{1\}\\leqDη​\(1−εNS\)​K\+LR​\(1\+εNS\)2​η1−εNS\+2​G1​ρ¯K\+LR​q¯Kη​\(1−εNS\)\\displaystyle\\frac\{D\}\{\\eta\(1\-\\varepsilon\_\{\\mathrm\{NS\}\}\)K\}\+\\frac\{L\_\{R\}\(1\+\\varepsilon\_\{\\mathrm\{NS\}\}\)^\{2\}\\eta\}\{1\-\\varepsilon\_\{\\mathrm\{NS\}\}\}\+\\frac\{2G\_\{1\}\\bar\{\\rho\}\_\{K\}\+L\_\{R\}\\bar\{q\}\_\{K\}\}\{\\eta\(1\-\\varepsilon\_\{\\mathrm\{NS\}\}\)\}\+\(1\+1\+εNS1−εNS\)​1K​∑k=0K−1𝔼​Tk\.\\displaystyle\+\\left\(1\+\\frac\{1\+\\varepsilon\_\{\\mathrm\{NS\}\}\}\{1\-\\varepsilon\_\{\\mathrm\{NS\}\}\}\\right\)\\frac\{1\}\{K\}\\sum\_\{k=0\}^\{K\-1\}\\mathbb\{E\}T\_\{k\}\.Lemma[5](https://arxiv.org/html/2608.06554#Thmlemma5), together with‖M~k‖∞≤\(1\+εNS\)\\\|\\widetilde\{M\}\_\{k\}\\\|\_\{\\infty\}\\leq\(1\+\\varepsilon\_\{\\mathrm\{NS\}\}\), gives

1K​∑k=0K−1𝔼​Tk≤𝔼​‖M1−H0‖1\(1−β\)​K\+n​σb\+β​LR1−β​\(η​\(1\+εNS\)\+ρ¯K\)\.\\frac\{1\}\{K\}\\sum\_\{k=0\}^\{K\-1\}\\mathbb\{E\}T\_\{k\}\\leq\\frac\{\\mathbb\{E\}\\\|M\_\{1\}\-H\_\{0\}\\\|\_\{1\}\}\{\(1\-\\beta\)K\}\+\\frac\{\\sqrt\{n\}\\,\\sigma\}\{\\sqrt\{b\}\}\+\\frac\{\\beta L\_\{R\}\}\{1\-\\beta\}\\left\(\\eta\(1\+\\varepsilon\_\{\\mathrm\{NS\}\}\)\+\\bar\{\\rho\}\_\{K\}\\right\)\.Substituting this estimate into the previous display proves the theorem\. ∎

## Acknowledgments

Portions of this work were conducted using the advanced computing resources provided by Texas A&M High Performance Research Computing \(HPRC\)\.

## References

- Adhikary et al\., \(2020\)Adhikary, S\., Srinivasan, S\., Gordon, G\., and Boots, B\. \(2020\)\.Expressiveness and learning of hidden quantum Markov models\.InProceedings of the Twenty Third International Conference on Artificial Intelligence and Statistics, volume 108 ofProceedings of Machine Learning Research, pages 4151–4161\. PMLR\.
- Agarwal et al\., \(2019\)Agarwal, V\., Reddy, N\. J\. K\., and Anand, A\. \(2019\)\.Unsupervised representation learning of DNA sequences\.arXiv preprint arXiv:1906\.03087\.
- Beck, \(2017\)Beck, A\. \(2017\)\.First\-order methods in optimization\.SIAM\.
- Biamonte et al\., \(2016\)Biamonte, J\., Wittek, P\., Pancotti, N\., Rebentrost, P\., Wiebe, N\., and Lloyd, S\. \(2016\)\.Quantum machine learning\.arXiv preprint arXiv:1611\.09347\.
- Burge and Karlin, \(1997\)Burge, C\. B\. and Karlin, S\. \(1997\)\.Modeling dependencies in pre\-mRNA splicing signals\.Proceedings of the National Academy of Sciences, 94\(2\):387–392\.
- Clark et al\., \(2015\)Clark, L\. A\., Huang, W\., Barlow, T\. M\., and Beige, A\. \(2015\)\.Hidden quantum Markov models and open quantum systems with instantaneous feedback\.InISCS 2014: Interdisciplinary Symposium on Complex Systems, pages 143–151\. Springer\.
- Deb et al\., \(2026\)Deb, B\., Pandey, S\., and Pakray, P\. \(2026\)\.Quantum hidden Markov model for sequential analysis\.Evolutionary Intelligence, 19:90\.
- Dheeru and Karra Taniskidou, \(2017\)Dheeru, D\. and Karra Taniskidou, E\. \(2017\)\.UCI machine learning repository\.[http://archive\.ics\.uci\.edu/ml](http://archive.ics.uci.edu/ml)\.
- Eddy, \(1998\)Eddy, S\. R\. \(1998\)\.Profile hidden Markov models\.Bioinformatics, 14\(9\):755–763\.
- Ghamsari, \(2025\)Ghamsari, M\. S\. \(2025\)\.Quantum computing applications in biology\.Discover Computing, 28\(1\):1–33\.
- Higham, \(1986\)Higham, N\. J\. \(1986\)\.Computing the polar decomposition—with applications\.SIAM Journal on Scientific and Statistical Computing, 7\(4\):1160–1174\.
- Higham, \(2008\)Higham, N\. J\. \(2008\)\.Functions of Matrices: Theory and Computation\.SIAM\.
- Jaeger, \(2000\)Jaeger, H\. \(2000\)\.Observable operator models for discrete stochastic time series\.Neural Computation, 12\(6\):1371–1398\.
- Jaggi, \(2013\)Jaggi, M\. \(2013\)\.Revisiting Frank\-Wolfe: Projection\-free sparse convex optimization\.InInternational conference on machine learning, pages 427–435\. PMLR\.
- Kim et al\., \(2026\)Kim, S\. W\. P\., von Keyserlingk, C\., and Lamacraft, A\. \(2026\)\.Measurement\-induced phase transitions in quantum inference problems and quantum hidden Markov models\.Physical Review Research, 8\(2\):023155\.
- Kösoglu\-Kind et al\., \(2023\)Kösoglu\-Kind, B\. et al\. \(2023\)\.A biological sequence comparison algorithm using quantum computers\.Scientific Reports, 13\(1\):14077\.
- Krogh et al\., \(1994\)Krogh, A\., Brown, M\., Mian, I\. S\., Sjolander, K\., and Haussler, D\. \(1994\)\.Hidden Markov models in computational biology: Applications to protein modeling\.Journal of Molecular Biology, 235\(5\):1501–1531\.
- Leifer and Poulin, \(2008\)Leifer, M\. S\. and Poulin, D\. \(2008\)\.Quantum graphical models and belief propagation\.Annals of Physics, 323\(8\):1899–1946\.
- Leifer and Spekkens, \(2013\)Leifer, M\. S\. and Spekkens, R\. W\. \(2013\)\.Towards a formulation of quantum theory as a causally neutral theory of Bayesian inference\.Physical Review A, 88\(5\):052130\.
- Monras et al\., \(2010\)Monras, A\., Beige, A\., and Wiesner, K\. \(2010\)\.Hidden quantum Markov models and non\-adaptive read\-out of many\-body states\.arXiv preprint arXiv:1002\.2337\.
- Nakatsukasa et al\., \(2010\)Nakatsukasa, Y\., Bai, Z\., and Gygi, F\. \(2010\)\.Optimizing Halley’s iteration for computing the matrix polar decomposition\.SIAM Journal on Matrix Analysis and Applications, 31\(5\):2700–2720\.
- Nalecz\-Charkiewicz et al\., \(2024\)Nalecz\-Charkiewicz, K\., Charkiewicz, K\., and Nowak, R\. M\. \(2024\)\.Quantum computing in bioinformatics: a systematic review mapping\.Briefings in Bioinformatics, 25\(5\):bbae391\.
- Nesterov, \(2013\)Nesterov, Y\. \(2013\)\.Introductory lectures on convex optimization: A basic course, volume 87\.Springer Science & Business Media\.
- Ning, \(2025\)Ning, N\. \(2025\)\.Robust iterative learning hidden quantum Markov models\.https://arxiv\.org/abs/2510\.23237\.
- Outeiral et al\., \(2021\)Outeiral, C\. et al\. \(2021\)\.The prospects of quantum computing in computational molecular biology\.WIREs Computational Molecular Science, 11\(1\):e1481\.
- Poddar et al\., \(2023\)Poddar, A\., Uzun, V\., Tunbridge, E\., Haerty, W\., and Nevado\-Holgado, A\. \(2023\)\.Identifying DNA sequence motifs using deep learning\.arXiv preprint arXiv:2311\.12884\.
- Schuld et al\., \(2015\)Schuld, M\., Sinayskiy, I\., and Petruccione, F\. \(2015\)\.An introduction to quantum machine learning\.Contemporary Physics, 56\(2\):172–185\.
- Souissi and Andolsi, \(2026\)Souissi, A\. and Andolsi, A\. \(2026\)\.A hidden quantum Markov model framework for entanglement and topological order in the AKLT chain\.The European Physical Journal Plus, 141\(4\):404\.
- \(29\)Srinivasan, S\., Downey, C\., and Boots, B\. \(2018a\)\.Learning and Inference in Hilbert space with Quantum Graphical Models\.InAdvances in Neural Information Processing Systems 31\.
- \(30\)Srinivasan, S\., Gordon, G\., and Boots, B\. \(2018b\)\.Learning hidden quantum Markov models\.InInternational Conference on Artificial Intelligence and Statistics, pages 1979–1987\.
- Towell and Shavlik, \(1991\)Towell, G\. G\. and Shavlik, J\. W\. \(1991\)\.Refinement of approximate domain theories by knowledge\-based neural networks\.Proceedings of the Eighth International Workshop on Machine Learning, pages 861–865\.
- Varsamis et al\., \(2023\)Varsamis, G\. D\. et al\. \(2023\)\.Quantum gate algorithm for reference\-guided DNA sequence alignment\.Briefings in Functional Genomics\.
- Warmuth and Kuzmin, \(2006\)Warmuth, M\. K\. and Kuzmin, D\. \(2006\)\.A bayesian probability calculus for density matrices\.InUAI\.
- Yeang, \(2010\)Yeang, C\.\-H\. \(2010\)\.A probabilistic graphical model of quantum systems\.In2010 Ninth International Conference on Machine Learning and Applications, pages 155–162\. IEEE\.
- Zhu et al\., \(2025\)Zhu, Q\.\-S\., Chen, G\., Yu, L\.\-H\., Xing, X\., and Li, X\.\-Y\. \(2025\)\.Channel\-constrained Markovian quantum diffusion model from open system perspective\.arXiv preprint arXiv:2511\.12221\.
- Åke Björck and Bowie, \(1971\)Åke Björck and Bowie, C\. \(1971\)\.An iterative algorithm for computing the best estimate of an orthogonal matrix\.SIAM Journal on Numerical Analysis, 8\(2\):358–364\.