LinOSS Stability

aThe referenced definitions, propositions, and corollaries are based on Jungers, d-LinOSS, and Bof.

LTI (Non-selective)¶

For LTI SSMs the transition matrix $M$ is a single fixed matrix across all inputs, so the maximal growth rate of the state $||x_{t}||$ is given by the spectral radius, i.e. the largest eigenvalue modulus (Gelfand's formula, Jungers Eq. 1.2):

\[ \rho(M)=\lim\limits_{t\to\infty}||M^{t}||^{\frac{1}{t}}=\max\limits_{i}|\lambda_{i}| \]

Two types of stability are of interest for SSMs:

Asymptotic Stability: Without force, the state decays to zero, for LTI systems iff $\rho(M)<1$.
Marginal (Lyapunov) Stability: Without force, the state remains bounded, iff $\rho(M)\le1$ with semisimple unit-modulus eigenvalues (e.g. undamped LinOSS-IMEX with $|\lambda_{i}|=1$).

For LinOSS in all its non-selective variants, Lyapunov Stability is proven by locating the eigenvalues of the fixed transition matrix in the unit disk. Since all parameters are diagonal, the oscillators uncouple (d-LinOSS Prop. A.1.), so we only look at the $i$-th oscillator and set $\det(\lambda I -[M]_{i})\stackrel{!}{=}0$. For the IMEX discretization derived below (d-LinOSS Prop. 3.1):

\[ \lambda^{\text{IMEX}}_{i_{1,2}}= \frac {1 + \frac{\Delta t_{i}}{2}G_{i} - \frac{\Delta t_{i}^{2}}{2}A_{i}} {1+\Delta t_{i}G_{i}} \pm \frac {\frac{\Delta t_{i}}{2}\sqrt{(G_{i}-\Delta t_{i}A_{i})^{2}-4A_{i}}} {1+\Delta t_{i}G_{i}} \]

Complex-conjugate pairs are guaranteed by the (sufficient) stability criterion $(G_{i}-\Delta t_{i}A_{i})^{2}\le4A_{i}$, making the discriminant non-positive. The magnitude of the complex pairs is then given by

\[ |\lambda^{\text{IMEX}}_{i}|=\frac{1}{\sqrt{1+\Delta t_{i}G_{i}}}\le1 \]

This holds for any non-negative $G_{i}$ and $\Delta t_{i}$ under the stability criterion. The same computation for the IM discretization gives

\[ \lambda^{\text{IM}}_{i_{1,2}}= \frac {1+\frac{\Delta t_{i}}{2}G_{i}\pm\frac{\Delta t_{i}}{2}\sqrt{G_{i}^{2}-4A_{i}}} {1+\Delta t_{i}G_{i}+\Delta t_{i}^{2}A_{i}}, \qquad |\lambda^{\text{IM}}_{i}|=\frac{1}{\sqrt{1+\Delta t_{i}G_{i}+\Delta t_{i}^{2}A_{i}}} \]

which is Schur for any $A_{i}>0$, $G_{i}\ge0$ with no criterion attached (the magnitude formula holding in the complex regime $G_{i}^{2}\le4A_{i}$). Already in the LTI case the two schemes split: IM is unconditionally contractive but its numerical dissipation entangles decay with frequency, while IMEX is dissipation-free at $G=0$ ($|\lambda|=1$, symplectic) and controls decay through $G$ alone, however requires constraining $A$ and $G$ to satisfy the stability criterion. This split is also observed under selectivity.

LTV (Selective)¶

Joint Spectral Radius¶

Selectivity can be analyzed in terms of the Joint Spectral Radius (JSR). Making any of the state matrix $A$, the damping $G$, or the time-step $\Delta t$ selective means we make the transition matrix $M$ depend on the input at time $t$. Thus, LinOSS, and selective SSMs in general, are an example of a switched dynamical system:

\[ \begin{align} x_{t+1} &= M_{t}x_{t}& M_{t} &\in \mathcal{M} \\ x_{0} &\in \mathbb{R}^{n} \end{align} \]

where $\mathcal{M}\subset\mathbb{R}^{n\times n}$ is the bounded set of parametric transition matrices. Such a system is uniformly exponentially stable under arbitrary switching if and only if the joint spectral radius is smaller than 1, i.e. $\rho(\mathcal{M})<1$ (Jungers Definition 1.1 and Corollary 1.1).

In case of a selective diagonal first-order model like Mamba, the discretization $\bar{A}=\exp(\Delta t A)$ with real negative $A$ gives a family of normal, commuting and simultaneously diagonalizable matrices, so a fixed eigenbasis across all inputs. In this case, the JSR is given by the largest eigenvalue magnitude of any matrix in the set $\mathcal{M}$ (Jungers Proposition 2.2) and it is then simply shown that $\rho(\mathcal{M})=\rho_{1}(\mathcal{M})<1$ like in the LTI case — provided the parameter set is compact with $\sup|\lambda|<1$, e.g. $\Delta t$ bounded below.

For LinOSS, if we let parameters vary arbitrarily, $\mathcal{M}$ is neither normal nor simultaneously triangularizable, so we can't apply the same logic as in the first-order diagonal system. Furthermore, the set of matrix families with $\rho(\mathcal{M})\le1$ is not semi-algebraic in general (Jungers Theorem 2.5), so an exact closed-form characterization in the style of d-LinOSS Proposition A.2. cannot be expected for the switched family; we settle for a sufficient criterion. To prove stability, we can leverage two facts:

$\mathcal{M}$ can be reduced into its subblocks (uncoupled oscillators, d-LinOSS Proposition A.1.), thus we only have to consider $\rho(\mathcal{M})=\max_{i}\rho({[M]_{i}})$ (Jungers Proposition 1.5). Since $A$, $G_{n}$, $\Delta t_{n}$ are diagonal, every per-oscillator block computation below involves only commuting factors, which we manipulate as scalars.
A common Lyapunov Function can prove stability (Jungers Proposition 2.8).

Throughout, the stiffness $A>0$ is kept static while the damping $G_{n}=G(u_{n})>0$ and the time-step $\Delta t_{n}=\Delta t(u_{n})>0$ are selective. This restriction is required for finding common metric $P$ below, which references $A$, and thus for finding a common Lyapunov Function. Time-varying metrics are possible, but selective stiffness may introduce Parametric Resonance and thus further constrain the system.

Continuous Time System and Energy¶

For d-LinOSS, the system is governed by

\[ \ddot{y}(t)=-A y(t)-G\dot{y}(t)+Bu(t) \]

which can be represented as two coupled first-order systems:

\[ \begin{align} \dot{z}(t) &= -Ay(t)-Gz(t)+Bu(t) \\ \dot{y}(t) &= z(t) \end{align} \]

The LinOSS discretizes a system of forced harmonic oscillators, which, undamped and unforced, conserve the energy $E(z,y)=\frac{1}{2}z^{2}+\frac{1}{2}ay^{2}$ (per oscillator, block reduction above). Differentiating along unforced trajectories:

\[ \begin{align} \dot{E} &= z\dot{z}+ay\dot{y} \\ &= z(-ay-gz)+ayz \\ &= -gz^{2}\le0 \end{align} \]

We have a dissipative Hamiltonian, which suggests $E$ as common Lyapunov function candidate, if positive-definiteness holds under discretization. The energy derivative and Lyapunov stability under both the IMEX and IM discretization can now be derived.

Selective IM Discretization¶

Both state variables are treated implicitly; the scheme inherits and amplifies the continuous dissipation.

IM Discretized State Matrix¶

Applying IM discretization:

\[ z_{n}=z_{n-1}+\Delta t_{n}(-Ay_{n}-G_{n}z_{n}+B_{n}u_{n}),\quad y_{n}=y_{n-1}+\Delta t_{n}z_{n} \]

or in matrix form with $x_{n}=[z_{n},y_{n}]^{\mathsf{T}}$:

\[ \begin{align} x_{n} &= x_{n-1}+ \begin{bmatrix} -\Delta t_{n} G_{n} & -\Delta t_{n}A \\ \Delta t_{n} \mathbf{I} & 0 \\ \end{bmatrix} x_{n} + \begin{bmatrix} \Delta t_{n}B_{n} \\ 0 \end{bmatrix}u_{n} \\[1em] \iff \begin{bmatrix} \mathbf{I} + \Delta t_{n} G_{n} & \Delta t_{n}A \\ -\Delta t_{n} \mathbf{I} & \mathbf{I} \\ \end{bmatrix} x_{n} &= x_{n-1}+ \begin{bmatrix} \Delta t_{n}B_{n} \\ 0 \end{bmatrix}u_{n} \end{align} \]

Inverting the left-hand matrix gives the IM discretized state matrix such that $x_{n}=M_{n}x_{n-1}+F_{n}u_{n}$:

\[ \begin{align} M_{n}&= \begin{bmatrix} S_{n} & -\Delta{t}_{n}S_{n}A \\ \Delta{t}_{n}S_{n} & \mathbf{I}-\Delta t_{n}^{2}S_{n}A \end{bmatrix} \\[1em] &= S_{n}\begin{bmatrix} \mathbf{I} & -\Delta{t}_{n}A \\ \Delta{t}_{n} & \mathbf{I}+\Delta{t}_{n}G_{n} \end{bmatrix} \end{align} \]

and

\[ F_{n}= S_{n}\begin{bmatrix} \Delta t_{n}B_{n} \\ \Delta t_{n}^{2}B_{n} \end{bmatrix} \]

with Schur complement $S^{-1}_{n}=\mathbf{I}+\Delta t_{n}G_{n} + \Delta t_{n}^{2}A$ and identity $\mathbf{I}-\Delta t_{n}^{2}S_{n}A=S_{n}(\mathbf{I}+\Delta t_{n}G_{n})$.

Discrete Time Energy Using IM Discretization¶

We now consider the energy evolution $\Delta{E}_{n} = z_{n}^{\mathsf{T}}\Delta{z}_{n}+y_{n}^{\mathsf{T}}A\Delta{y}_{n}$ under the IM update rule with the new states $z_{n}$ and $Ay_{n}$ (variables are treated implicitly), using the polarization identity $2b^{\mathsf{T}}(b-c)=|b|^{2}-|c|^{2}+|b-c|^{2}$:

$$

\begin{align} z_{n}^{\mathsf{T}}(z_{n}-z_{n-1}) &= \tfrac{1}{2}\big(|z_{n}|^{{2}-|z_{n-1}|} \[1em]}\big)+\tfrac{1}{2}|\Delta z_{n}|^{2} = -\Delta t_{n}z_{n}^{\mathsf{T}}Ay_{n}-\Delta t_{n}z_{n}^{\mathsf{T}}G_{n}z_{n

y_{n}^{\mathsf{T}}A(y_{n}-y_{n-1}) &= \tfrac{1}{2}\big(|y_{n}|{A}^{{2}-|y_{n-1}|_{A}}|\Delta y}\big)+\tfrac{1}{2|{A}^{2} = \Delta t \end{align}}z_{n}^{\mathsf{T}}Ay_{n

$$

The stiffness cross-terms cancel as in continuous time, and both increment terms enter with a negative sign. Summing and substituting $\Delta y_{n}=\Delta t_{n}z_{n}$:

\[ \Delta E_{n}= -z_{n}^{\mathsf{T}}\Big(\Delta t_{n}G_{n}+\tfrac{\Delta t_{n}^{2}}{2}A\Big)z_{n}-\tfrac{\Delta t_{n}^{2}}{2}\big|Ay_{n}+G_{n}z_{n}\big|^{2}\le0 \]

This is the discrete counterpart of $\dot{E}=-gz^{2}$: the first term is the physical damping, the remainder is the numerical dissipation of the implicit method. $E$ decreases for any $G_{n}**\ge**0$ and $\Delta t_{n}>0$ — even undamped — so $E$ is a common Lyapunov candidate with no criterion attached.

Lyapunov Stability in IM Discretization¶

Following Bof et al., choose the quadratic Lyapunov Function $V=2E(x)=x^{\mathsf{T}}Px$ using

\[P=\begin{bmatrix} \mathbf{I} & 0 \\ 0 & A \end{bmatrix}\]

Since $A$ is static, $P$ is constant; it is symmetric and bounded positive for $0<a_{\min}\le A\le a_{\max}$. We check stability by computing

\[ Q_{n}=P-M_{n}^{\mathsf{T}}PM_{n} \]

and asserting $Q_{n}$ to be symmetric and positive definite.

$$

\begin{align} M_{n}^{\mathsf{T}}PM_{n}

&= S_{n}\begin{bmatrix} \mathbf{I} & \Delta t_{n} \ -\Delta t_{n}A & \mathbf{I}+\Delta t_{n}G_{n} \end{bmatrix} \begin{bmatrix} \mathbf{I} & 0 \ 0 & A \end{bmatrix} S_{n} \begin{bmatrix} \mathbf{I} & -\Delta t_{n}A \ \Delta t_{n} & \mathbf{I}+\Delta t_{n}G_{n} \end{bmatrix} \[1em]

&= S_{n}^{2}\begin{bmatrix} \mathbf{I} & \Delta t_{n}A \ -\Delta t_{n}A & A(\mathbf{I}+\Delta t_{n}G_{n}) \end{bmatrix} \begin{bmatrix} \mathbf{I} & -\Delta t_{n}A \ \Delta t_{n} & \mathbf{I}+\Delta t_{n}G_{n} \end{bmatrix} \[1em] &= S_{n}^{2}\begin{bmatrix} \mathbf{I}+\Delta t_{n}^{2}A & -\Delta t_{n}A+\Delta t_{n}A+\Delta t_{n}^{2}AG_{n} \ -\Delta t_{n}A+\Delta t_{n}A+\Delta t_{n}^{2}AG_{n} & \Delta t_{n}^{2}A}+A+2\Delta t_{n}AG_{n}+\Delta t_{n^{2}AG_{n} \end{bmatrix} \[1em] &= S_{n}^{2}\begin{bmatrix} \mathbf{I}+\Delta t_{n}^{2}A & \Delta t_{n}^{2}AG_{n} \ \Delta t_{n}^{2}AG_{n} & \Delta t_{n}^{2}A \end{bmatrix} \end{align} $$}+A(\mathbf{I}+ \Delta t_{n}G_{n})^{2

Then

$$ \begin{align} Q_{n} &= P - M_{n}^{\mathsf{T}}PM_{n} \[1em]

&= \begin{bmatrix} \mathbf{I} & 0 \ 0 & A \end{bmatrix} - S_{n}^{2}\begin{bmatrix} \mathbf{I}+\Delta t_{n}^{2}A & \Delta t_{n}^{2}AG_{n} \ \Delta t_{n}^{2}AG_{n} & \Delta t_{n}^{2}A \end{bmatrix} \[1em] &= \begin{bmatrix} \mathbf{I} - S_{n}^{2}(\mathbf{I} + \Delta t_{n}^{2}A) & -S_{n}^{2}\Delta t_{n}^{2}AG_{n} \ -S_{n}^{2}\Delta t_{n}^{2}AG_{n} & A( \mathbf{I} - S_{n}^{2}[ \Delta t_{n}^{2}A + (\mathbf{I} + \Delta t_{n}G_{n})^{2} ] ) \end{bmatrix} \end{align} $$}+A(\mathbf{I}+ \Delta t_{n}G_{n})^{2

Given $S_{n}^{-1} = \mathbf{I} + \Delta t_{n}G_{n} + \Delta t_{n}^{2}A$, we simplify the blocks of $Q_{n}$:

$$ \begin{aligned} Q_{11} &= \mathbf{I} - S_{n}^{2}(\mathbf{I} + \Delta t_{n}^{2}A) \ &= \mathbf{I} - S_{n}^{2}(S_{n}) \ &= \mathbf{I} - S_{n} + S_{n}^{2}\Delta t_{n}G_{n} \[1em]} - \Delta t_{n}G_{n

Q_{22} &= A ( \mathbf{I} - S_{n}^{2}[ \Delta t_{n}^{2}A + (\mathbf{I} + \Delta t_{n}G_{n})^{2} ] ) \ &= A ( \mathbf{I} - S_{n}^{2}[ S_{n}^{-1} + \Delta t_{n}G_{n}(\mathbf{I} + \Delta t_{n}G_{n}) ] ) \ &= A ( \mathbf{I} - S_{n} - S_{n}^{2}\Delta t_{n}G_{n}(\mathbf{I} + \Delta t_{n}G_{n}) ) \end{aligned} $$

Thus, the simplified matrix $Q_{n}$ is:

\[ \begin{aligned} Q_{n} &= \begin{bmatrix} \mathbf{I} - S_{n} + S_{n}^{2}\Delta t_{n}G_{n} & -S_{n}^{2}\Delta t_{n}^{2}AG_{n} \\ -S_{n}^{2}\Delta t_{n}^{2}AG_{n} & A( \mathbf{I} - S_{n} - S_{n}^{2}\Delta t_{n}G_{n}(\mathbf{I} + \Delta t_{n}G_{n}) ) \end{bmatrix} \end{aligned} \]

For the determinant, substituting $S_{n}^{-1}$ into the block expansion and using $(\mathbf{I}+\Delta t_{n}^{2}A)[\Delta t_{n}^{2}A+(\mathbf{I}+\Delta t_{n}G_{n})^{2}]-\Delta t_{n}^{4}AG_{n}^{2}=(S_{n}^{-1})^{2}$ yields:

\[ \begin{aligned} \det(Q_{n}) &= \det( A { \mathbf{I} - S_{n}^{2} [ \mathbf{I} + 2\Delta t_{n}^{2}A + (\mathbf{I} + \Delta t_{n}G_{n})^{2} ] + S_{n}^{4} (S_{n}^{-1})^{2} } ) \\[1em] &= \det( A { \mathbf{I} - S_{n}^{2} [ 2S_{n}^{-1} + \Delta t_{n}^{2}G_{n}^{2} ] + S_{n}^{2} } ) \\[1em] &= \det( A { (\mathbf{I} - S_{n})^{2} - (S_{n}\Delta t_{n}G_{n})^{2} } ) \\[1em] &= \det( A [ \mathbf{I} - S_{n}(\mathbf{I} + \Delta t_{n}G_{n}) ] [ \mathbf{I} - S_{n}(\mathbf{I} - \Delta t_{n}G_{n}) ] ) \\[1em] &= \det( A [ S_{n}\Delta t_{n}^{2}A ] [ S_{n}(2\Delta t_{n}G_{n} + \Delta t_{n}^{2}A) ] ) \\[1em] &= \det( S_{n}^{2}\Delta t_{n}^{2}A^{2} ( 2\Delta t_{n}G_{n} + \Delta t_{n}^{2}A ) ) \end{aligned} \]

Since all blocks are commuting positive diagonals, $Q_{n}>0$ iff $Q_{11}>0$ and $\det(Q_{n})>0$ per oscillator. The determinant is positive for any $G_{n}\ge0$ whenever $A>0$ and $\Delta t_{n}>0$, and $Q_{11}=S_{n}(\Delta t_{n}G_{n}+\Delta t_{n}^{2}A)+S_{n}^{2}\Delta t_{n}G_{n}>0$ follows. Hence, mirroring the energy identity:

\[ Q_{n}>0 \quad\text{for all}\quad G_{n}\ge0,\ \Delta t_{n}>0 \]

Two regimes:

Parameters in a compact set ($0<\delta_{\min}\le\Delta t_{n}\le\delta_{\max}$, $0<a_{\min}\le A\le a_{\max}$, $G_{n}\le g_{\max}\mathbf{I}$): then $Q_{n}\ge q_{1}\mathbf{I}>0$ uniformly and the system is globally exponentially stable (Bof Theorem 5.9); the common Lyapunov function bounds the JSR, $\rho(\mathcal{M})<1$ (Jungers Proposition 2.8).
If we allow degenerate $\Delta t_{n},A\ge0$, the determinant may vanish and $Q_{n}$ is only positive semidefinite: uniformly stable but not necessarily convergent to equilibrium (Bof Theorem 5.4).

With forcing, $\Delta V_{n}\le -x_{n-1}^{\mathsf{T}}Q_{n}x_{n-1}+2\Delta t_{n}z_{n}^{\mathsf{T}}B_{n}u_{n}$, so Young's inequality gives an ISS bound: bounded inputs yield bounded states.

So we have a static bank of oscillators with input-dependent time-evolution and damping, unconditionally stable: the implicit stiffness converts the resonance channel into numerical dissipation — at the cost of entangling decay with frequency, $|\lambda_{n}|=(1+\Delta t_{n}G_{n}+\Delta t_{n}^{2}A)^{-1/2}$.

Selective IMEX Discretization¶

The stiffness is treated explicitly (d-LinOSS Eq. (4)): IMEX trades IM's numerical dissipation for exact energy behavior — the undamped scheme is symplectic — and the dissipation balance becomes conditional.

IMEX Discretized State Matrix¶

Applying IMEX discretization with per-step parameters:

\[ z_{n}=z_{n-1}+\Delta t_{n}(-Ay_{n-1}-G_{n}z_{n}+B_{n}u_{n}),\quad y_{n}=y_{n-1}+\Delta t_{n}z_{n} \]

or in matrix form:

$$ \begin{bmatrix} \mathbf{I}+\Delta t_{n}G_{n} & 0 \ -\Delta t_{n}\mathbf{I} & \mathbf{I} \end{bmatrix} x_{n}

= \begin{bmatrix} \mathbf{I} & -\Delta t_{n}A \ 0 & \mathbf{I} \end{bmatrix} x_{n-1} +\begin{bmatrix} \Delta t_{n}B_{n} \ 0 \end{bmatrix}u_{n} $$

The left-hand matrix is block lower-triangular, so only the Schur complement $S_{n}^{-1}=\mathbf{I}+\Delta t_{n}G_{n}$ must be inverted (contrast IM: $S^{-1}=\mathbf{I}+\Delta t G+\Delta t^{2}A$). Inverting gives $x_{n}=M_{n}x_{n-1}+F_{n}u_{n}$ with

\[ M_{n}=\begin{bmatrix} S_{n} & -\Delta t_{n}S_{n}A \\ \Delta t_{n}S_{n} & \mathbf{I}-\Delta t_{n}^{2}S_{n}A \end{bmatrix} = S_{n}\begin{bmatrix} \mathbf{I} & -\Delta t_{n}A \\ \Delta t_{n} & S_{n}^{-1}-\Delta t_{n}^{2}A \end{bmatrix} \]

and

\[ F_{n}= S_{n}\begin{bmatrix} \Delta t_{n}B_{n} \\ \Delta t_{n}^{2}B_{n} \end{bmatrix} \]

using the identity $\mathbf{I}-\Delta t_{n}^{2}S_{n}A=S_{n}(S_{n}^{-1}-\Delta t_{n}^{2}A)$. This is d-LinOSS Eq. (6) with time indices, and the frozen-step spectrum is the LTI IMEX spectrum above with time indices: under $(G_{n}-\Delta t_{n}A)^{2}\le 4A$, $|\lambda_{n}|=(1+\Delta t_{n}G_{n})^{-1/2}\le1$. However, since the $M_{n}$ do not commute, per-step Schur stability certifies nothing about $\rho(\mathcal{M})$ — see the failure mode below.

Discrete Time Energy Using IMEX Discretization¶

Test the $z$-update with the new state $z_{n}$ and the $y$-update with the old state $Ay_{n-1}$ (the explicitly treated variable), using the polarization identity $2b^{\mathsf{T}}(b-c)=|b|^{2}-|c|^{2}+|b-c|^{2}$:

\[ \begin{align} z_{n}^{\mathsf{T}}(z_{n}-z_{n-1}) &= \tfrac{1}{2}\big(|z_{n}|^{2}-|z_{n-1}|^{2}\big)+\tfrac{1}{2}|\Delta z_{n}|^{2} = -\Delta t_{n}z_{n}^{\mathsf{T}}Ay_{n-1}-\Delta t_{n}z_{n}^{\mathsf{T}}G_{n}z_{n} \\ y_{n-1}^{\mathsf{T}}A(y_{n}-y_{n-1}) &= \tfrac{1}{2}\big(|y_{n}|_{A}^{2}-|y_{n-1}|_{A}^{2}\big)-\tfrac{1}{2}|\Delta y_{n}|_{A}^{2} = \Delta t_{n}z_{n}^{\mathsf{T}}Ay_{n-1} \end{align} \]

The stiffness cross-terms cancel as in continuous time, but testing the explicit variable at the old state flips the sign of its increment term. Summing and substituting $\Delta y_{n}=\Delta t_{n}z_{n}$:

\[ \Delta E_{n}= -z_{n}^{\mathsf{T}}\Big(\Delta t_{n}G_{n}-\tfrac{\Delta t_{n}^{2}}{2}A\Big)z_{n}-\tfrac{\Delta t_{n}^{2}}{2}\big|Ay_{n-1}+G_{n}z_{n}\big|^{2} \]

Unlike IM, the explicit stiffness injects energy $+\tfrac{\Delta t_{n}^{2}}{2}z_{n}^{\mathsf{T}}Az_{n}$ per step, which the damping must absorb. This suggests $E$ remains a Lyapunov candidate only above a damping floor $G_{n}\ge\tfrac{\Delta t_{n}}{2}A$, made precise next.

Lyapunov Stability in IMEX Discretization¶

Following Bof et al., choose again the quadratic Lyapunov Function $V=2E(x)=x^{\mathsf{T}}Px$ with the same constant metric $P$. We check stability by computing

\[ Q_{n}=P-M_{n}^{\mathsf{T}}PM_{n} \]

and asserting $Q_{n}$ to be symmetric and positive definite.

\[ \begin{align} M_{n}^{\mathsf{T}}PM_{n} &= S_{n}\begin{bmatrix} \mathbf{I} & \Delta t_{n} \\ -\Delta t_{n}A & S_{n}^{-1}-\Delta t_{n}^{2}A \end{bmatrix} \begin{bmatrix} \mathbf{I} & 0 \\ 0 & A \end{bmatrix} S_{n}\begin{bmatrix} \mathbf{I} & -\Delta t_{n}A \\ \Delta t_{n} & S_{n}^{-1}-\Delta t_{n}^{2}A \end{bmatrix} \\[1em] &= S_{n}^{2}\begin{bmatrix} \mathbf{I} & \Delta t_{n}A \\ -\Delta t_{n}A & A,(S_{n}^{-1}-\Delta t_{n}^{2}A) \end{bmatrix} \begin{bmatrix} \mathbf{I} & -\Delta t_{n}A \\ \Delta t_{n} & S_{n}^{-1}-\Delta t_{n}^{2}A \end{bmatrix} \\[1em] &= S_{n}^{2}\begin{bmatrix} \mathbf{I}+\Delta t_{n}^{2}A & -\Delta t_{n}A+\Delta t_{n}A(S_{n}^{-1}-\Delta t_{n}^{2}A) \\ -\Delta t_{n}A+\Delta t_{n}A(S_{n}^{-1}-\Delta t_{n}^{2}A) & \Delta t_{n}^{2}A^{2}+A(S_{n}^{-1}-\Delta t_{n}^{2}A)^{2} \end{bmatrix} \\[1em] &= S_{n}^{2}\begin{bmatrix} \mathbf{I}+\Delta t_{n}^{2}A & \Delta t_{n}^{2}A(G_{n}-\Delta t_{n}A) \\ \Delta t_{n}^{2}A(G_{n}-\Delta t_{n}A) & \Delta t_{n}^{2}A^{2}+A(S_{n}^{-1}-\Delta t_{n}^{2}A)^{2} \end{bmatrix} \end{align} \]

using $S_{n}^{-1}-\mathbf{I}=\Delta t_{n}G_{n}$ in the off-diagonal blocks. Then

\[ Q_{n}= \begin{bmatrix} \mathbf{I}-S_{n}^{2}(\mathbf{I}+\Delta t_{n}^{2}A) & -S_{n}^{2}\Delta t_{n}^{2}A(G_{n}-\Delta t_{n}A) \\ -S_{n}^{2}\Delta t_{n}^{2}A(G_{n}-\Delta t_{n}A) & A(\mathbf{I}-S_{n}^{2}[\Delta t_{n}^{2}A+(S_{n}^{-1}-\Delta t_{n}^{2}A)^{2}]) \end{bmatrix} \]

Given $S_{n}^{-1}=\mathbf{I}+\Delta t_{n}G_{n}$, we simplify the blocks of $Q_{n}$:

\[ \begin{aligned} Q_{11} &= \mathbf{I}-S_{n}^{2}(\mathbf{I}+\Delta t_{n}^{2}A) = S_{n}^{2}\big(S_{n}^{-2}-\mathbf{I}-\Delta t_{n}^{2}A\big) \\ &= S_{n}^{2}\big((\mathbf{I}+\Delta t_{n}G_{n})^{2}-\mathbf{I}-\Delta t_{n}^{2}A\big) \\ &= S_{n}^{2},\Delta t_{n}\big(2G_{n}+\Delta t_{n}(G_{n}^{2}-A)\big) \\[1em] Q_{22} &= A(\mathbf{I}-S_{n}^{2}[\Delta t_{n}^{2}A+S_{n}^{-2}-2\Delta t_{n}^{2}S_{n}^{-1}A+\Delta t_{n}^{4}A^{2}]) \\ &= A,\Delta t_{n}^{2}S_{n}A\big(2\mathbf{I}-S_{n}(\mathbf{I}+\Delta t_{n}^{2}A)\big) \\ &= S_{n}^{2},\Delta t_{n}^{2}A^{2}\big(\mathbf{I}+2\Delta t_{n}G_{n}-\Delta t_{n}^{2}A\big) \end{aligned} \]

where the last line uses $2S_{n}^{-1}-\mathbf{I}-\Delta t_{n}^{2}A=\mathbf{I}+2\Delta t_{n}G_{n}-\Delta t_{n}^{2}A$. Thus, the simplified matrix $Q_{n}$ is:

\[ Q_{n}=S_{n}^{2} \begin{bmatrix} \Delta t_{n}\big(2G_{n}+\Delta t_{n}(G_{n}^{2}-A)\big) & -\Delta t_{n}^{2}A(G_{n}-\Delta t_{n}A) \\ -\Delta t_{n}^{2}A(G_{n}-\Delta t_{n}A) & \Delta t_{n}^{2}A^{2}\big(\mathbf{I}+2\Delta t_{n}G_{n}-\Delta t_{n}^{2}A\big) \end{bmatrix} \]

For the determinant, expanding the block products and collecting terms yields:

\[ \begin{aligned} \det(Q_{n}) &= \det\Big( S_{n}^{4}\Delta t_{n}^{3}A^{2}\Big[ \big(2G_{n}+\Delta t_{n}G_{n}^{2}-\Delta t_{n}A\big)\big(\mathbf{I}+2\Delta t_{n}G_{n}-\Delta t_{n}^{2}A\big)-\Delta t_{n}\big(G_{n}-\Delta t_{n}A\big)^{2} \Big] \Big) \\[1em] &= \det\Big( S_{n}^{4}\Delta t_{n}^{3}A^{2},\big(2G_{n}-\Delta t_{n}A\big)\big(\mathbf{I}+\Delta t_{n}G_{n}\big)^{2} \Big) \\[1em] &= \det\Big( S_{n}^{2},\Delta t_{n}^{3}A^{2},\big(2G_{n}-\Delta t_{n}A\big) \Big) \end{aligned} \]

since the bracket factorizes as $(2G_{n}-\Delta t_{n}A)(\mathbf{I}+2\Delta t_{n}G_{n}+\Delta t_{n}^{2}G_{n}^{2})$ and $S_{n}^{4}(\mathbf{I}+\Delta t_{n}G_{n})^{2}=S_{n}^{2}$. Compare IM: $\det(Q_{n})=\det\big(S_{n}^{2}\Delta t_{n}^{2}A^{2}(2\Delta t_{n}G_{n}+\Delta t_{n}^{2}A)\big)$, which is positive for any $G_{n}\ge0$ — here the factor $(2G_{n}-\Delta t_{n}A)$ can change sign.

Since all blocks are commuting positive diagonals, $Q_{n}>0$ iff $Q_{11}>0$ and $\det(Q_{n})>0$ per oscillator. The determinant is positive iff $2G_{n}-\Delta t_{n}A>0$, and under this condition $Q_{11}=S_{n}^{2}\Delta t_{n}\big(\Delta t_{n}G_{n}^{2}+(2G_{n}-\Delta t_{n}A)\big)>0$ follows automatically. Hence, mirroring the energy identity:

\[ Q_{n}>0 \iff G_{n}>\tfrac{\Delta t_{n}}{2}A \]

i.e. per oscillator the damping floor $g_{i}(u_{n})>\tfrac{\Delta t_{i}(u_{n})}{2}a_{i}$. Three regimes:

$G_{n}\ge\tfrac{\Delta t_{n}}{2}A+\epsilon\mathbf{I}$ with parameters in a compact set (as for IM): then $Q_{n}\ge q_{1}\mathbf{I}>0$ uniformly and the system is globally exponentially stable (Bof Theorem 5.9); the common Lyapunov function bounds the JSR, $\rho(\mathcal{M})<1$ (Jungers Proposition 2.8).
$G_{n}=\tfrac{\Delta t_{n}}{2}A$ exactly: $\det(Q_{n})=0$ and $Q_{11}=S_{n}^{2}\Delta t_{n}^{2}G_{n}^{2}\ge0$, so $Q_{n}\ge0$: uniformly stable but not necessarily convergent (Bof Theorem 5.4).
$G_{n}<\tfrac{\Delta t_{n}}{2}A$ somewhere: $Q_{n}$ indefinite, even if $M_{n}$ is Schur (e.g. $A=1$, $\Delta t=0.8$, $G=0.3$: $|\lambda|=0.898$ but $\lambda(Q)={-0.15,,0.44}$) — and instability under switching is real, not an artifact of the choice of $P$ (next section).

Enforcing the damping floor is enough for stability, thus the d-LinOSS stability criterion is optional: - Keeping it restricts eigenvalues to unit-disk and gives a bijection between parameter space and eigenvalues in the unit-disk. - Removing it lets eigenvalues become real, enabling overdamped system dynamics (real decay at variable rates).

With forcing, $\Delta V_{n}\le -x_{n-1}^{\mathsf{T}}Q_{n}x_{n-1}+2\Delta t_{n}z_{n}^{\mathsf{T}}B_{n}u_{n}$, so above the floor Young's inequality gives an ISS bound: bounded inputs yield bounded states.

Failure Below the Floor: Parametric Resonance¶

Per-step spectral stability does not compose. With $A=1$, $G_{n}=0$ and alternating $\Delta t_{n}\in\{0.25,0.95\}$, each matrix satisfies the d-LinOSS criterion with $|\lambda_{n}|=1$, yet

\[ \rho(M_{2}M_{1})^{1/2}\approx1.565 \]

per step (worst alternating pair on $(0,2]^{2}$: $\approx1.72$); even i.i.d. random $\Delta t_{n}$ grows at $\approx1.07$ per step. Mechanism: each undamped IMEX step is symplectic ($\det[M_{n}]_{i}=1$) and preserves its own $\Delta t$-dependent modified Hamiltonian; switching $\Delta t$ switches the invariant ellipse and pumps energy — discrete Parametric Resonance through the time-step alone. So $\rho(\mathcal{M})>1$ for the unconstrained family, and the floor is nearly sharp: greedy adversarial switching still finds growth at half the floor.

So contrary to IM (unconditionally stable, at the cost of $|\lambda|=(1+\Delta t G+\Delta t^{2}A)^{-1/2}$ entangling decay with frequency), the selective IMEX bank of oscillators with input-dependent damping and time-evolution is stable iff the damping absorbs the explicit-stiffness energy injection, $G_{n}\ge\tfrac{\Delta t_{n}}{2}A$, which caps

\[ |\lambda_{n}|\le\Big(1+\tfrac{\Delta t_{n}^{2}A}{2}\Big)^{-1/2} \]

which is half the numerical dissipation of IM's floor. The constraint is enforceable by reparameterization in the style of d-LinOSS Eq. (8): $G_{n}=\tfrac{\Delta t_{n}}{2}A+\operatorname{softplus}(\cdot)$.