diff --git a/UNKNOWN.egg-info/SOURCES.txt b/UNKNOWN.egg-info/SOURCES.txt
new file mode 100644
index 0000000..0030353
--- /dev/null
+++ b/UNKNOWN.egg-info/SOURCES.txt
@@ -0,0 +1,7 @@
+LICENSE
+README.md
+pyproject.toml
+UNKNOWN.egg-info/PKG-INFO
+UNKNOWN.egg-info/SOURCES.txt
+UNKNOWN.egg-info/dependency_links.txt
+UNKNOWN.egg-info/top_level.txt
\ No newline at end of file
diff --git a/modpods/__init__.py b/modpods/__init__.py
index 60cb837..f1c81a2 100644
--- a/modpods/__init__.py
+++ b/modpods/__init__.py
@@ -3,8 +3,12 @@
 from .estimator import DelayIO, DelayIOModel
 from .kernels import (
     BimodalGammaKernel,
+    CanonicalLTIKernel,
     ConvolutionKernel,
+    DirectLTISystem,
+    ExponentialDecayKernel,
     ExponentialGrowthKernel,
+    ExponentialKernel,
     GammaKernel,
     LogNormalKernel,
     UnderdampedOscillatorKernel,
@@ -25,7 +29,7 @@
 from .model import SINDY_delays_MI
 from .predict import delay_io_predict
 from .topology import TopologyInference, find_topology_no_geo, infer_causative_topology
-from .train import delay_io_train
+from .train import decoupled_lti_train, delay_io_train
 from .transforms import (
     TransformCache,
     make_kernel_params,
@@ -40,10 +44,14 @@
     "DelayIO",
     "DelayIOModel",
     "ConvolutionKernel",
+    "CanonicalLTIKernel",
+    "DirectLTISystem",
     "GammaKernel",
     "LogNormalKernel",
     "BimodalGammaKernel",
+    "ExponentialDecayKernel",
     "ExponentialGrowthKernel",
+    "ExponentialKernel",
     "UnderdampedOscillatorKernel",
     "get_kernel",
     "list_kernels",
@@ -53,6 +61,8 @@
     "params_vector_to_dataframe",
     "transform_inputs",
     "delay_io_train",
+    "decoupled_lti_train",
+    "direct_lti_train",
     "SINDY_delays_MI",
     "delay_io_predict",
     "lti_from_gamma",
diff --git a/modpods/_system_id.py b/modpods/_system_id.py
index f764487..9970870 100644
--- a/modpods/_system_id.py
+++ b/modpods/_system_id.py
@@ -253,6 +253,7 @@ def _active_set_qp(
     d: np.ndarray,
     max_iter: int = 50,
     tol: float = 1e-8,
+    ridge_lambda: float = 1e-8,
 ) -> np.ndarray:
     """Solve  min ||A w - b||^2  s.t.  C w <= d  via the active-set method.
 
@@ -261,7 +262,10 @@ def _active_set_qp(
     available — cvxpy is an explicit dependency already.
     """
     n = A.shape[1]
-    w = np.linalg.lstsq(A, b, rcond=None)[0]
+    # Use regularized least squares for better numerical stability
+    AtA = A.T @ A + ridge_lambda * np.eye(n)
+    Atb = A.T @ b
+    w = np.linalg.solve(AtA, Atb)
     active: set[int] = set()
 
     for _ in range(max_iter):
@@ -277,10 +281,10 @@ def _active_set_qp(
         d_active = d[list(active)]
 
         # Equality-constrained least-squares via Lagrange multipliers
-        AtA = A.T @ A + np.eye(n) * 1e-10
-        Atb = A.T @ b
-        w_ls = np.linalg.solve(AtA, Atb)
-        A_inv = np.linalg.inv(AtA)
+        AtA_reg = A.T @ A + ridge_lambda * np.eye(n)
+        Atb_reg = A.T @ b
+        w_ls = np.linalg.solve(AtA_reg, Atb_reg)
+        A_inv = np.linalg.inv(AtA_reg)
         CAt = C_active @ A_inv
         denom = CAt @ C_active.T
         if denom.size == 1:
@@ -467,7 +471,7 @@ def fit(
         theta_valid = theta[valid]
         x_dot_valid = x_dot_arr[valid]
 
-        # Solve
+        # Solve with regularization
         self._coef = self._solve(theta_valid, x_dot_valid)
 
         # Cache computed arrays for potential reuse in score()
@@ -483,8 +487,12 @@ def fit(
     def _solve(self, theta: np.ndarray, x_dot: np.ndarray) -> np.ndarray:
         """Return coefficient matrix of shape (n_targets, n_features)."""
         if self.constraint_lhs is None or self.constraint_rhs is None:
-            # Unconstrained OLS (equivalent to STLSQ threshold=0, alpha=0)
-            coef = np.asarray(np.linalg.lstsq(theta, x_dot, rcond=None)[0])
+            # Regularized OLS (ridge regression) for better numerical stability
+            # This avoids SVD convergence issues with ill-conditioned matrices
+            ridge_lambda = 1e-8
+            AtA = theta.T @ theta + ridge_lambda * np.eye(theta.shape[1])
+            Atb = theta.T @ x_dot
+            coef = np.linalg.solve(AtA, Atb)
             return coef.T
         else:
             C = self.constraint_lhs
@@ -506,7 +514,9 @@ def _solve_equality_constrained(
         n_targets = x_dot.shape[1] if x_dot.ndim > 1 else 1
         x_dot_2d = x_dot.reshape(-1, 1) if x_dot.ndim == 1 else x_dot
 
-        AtA = theta.T @ theta + np.eye(n_feat) * 1e-10
+        # Add regularization for numerical stability
+        ridge_lambda = 1e-8
+        AtA = theta.T @ theta + ridge_lambda * np.eye(n_feat)
         Atb = theta.T @ x_dot_2d  # (n_feat, n_targets)
         w_ls = np.linalg.solve(AtA, Atb)  # (n_feat, n_targets)
         A_inv = np.linalg.inv(AtA)
diff --git a/modpods/kernels.py b/modpods/kernels.py
index dd300e8..c081663 100644
--- a/modpods/kernels.py
+++ b/modpods/kernels.py
@@ -12,6 +12,7 @@
 
 import numpy as np
 import scipy.stats as stats
+from scipy.linalg import expm
 
 
 class ConvolutionKernel(ABC):
@@ -65,8 +66,21 @@ def kernel_fn(self, t: np.ndarray, *params: float) -> np.ndarray:
         """
         ...
 
+    @property
+    def is_unstable(self) -> bool:
+        """Whether this kernel represents an unstable impulse response."""
+        return False
+
+    def is_unstable_params(self, *params: float) -> bool:
+        return self.is_unstable
+
+    def is_stable_delay(self, *params: float) -> bool:
+        return True
+
+    def to_lti(self, *params: float) -> tuple:
+        return None
+
     def make_kwargs(self, params: np.ndarray) -> dict:
-        """Convert flat parameter array to a kwargs dict keyed by param_names."""
         return dict(zip(self.param_names, params.tolist()))
 
 
@@ -201,7 +215,7 @@ class UnderdampedOscillatorKernel(ConvolutionKernel):
 
     Parameters are physical: zeta (damping ratio) and omega_n (natural frequency).
     Positive zeta produces decaying oscillations; negative zeta produces growing
-    (unstable) oscillations.  The kernel is truncated to non-negative values for
+    (unstable) oscillations. The kernel is truncated to non-negative values for
     causality when zeta >= 0.
 
     Note: This does NOT construct LTI state-space matrices.  It only uses the
@@ -225,8 +239,8 @@ def param_names(self) -> List[str]:
     def default_bounds(self) -> np.ndarray:
         return np.array(
             [
-                [-0.99, 0.99],
-                [0.1, 10.0],
+                [0.001, 5.0],
+                [0.001, 50.0],
             ]
         )
 
@@ -235,18 +249,44 @@ def default_init(self) -> np.ndarray:
         return np.array([0.1, 2.0])
 
     def kernel_fn(self, t: np.ndarray, zeta: float, omega_n: float) -> np.ndarray:  # type: ignore[override]
-        omega_d = omega_n * np.sqrt(1.0 - zeta**2)
-        amplitude = omega_n / omega_d
-        # For numerical stability, clip the exponent
-        exponent = -zeta * omega_n * t
-        # Clip exponent to prevent overflow (exp(700) ~ 1e304, near float64 max)
-        max_exponent = 700.0
-        exponent = np.clip(exponent, -max_exponent, max_exponent)
-        h = amplitude * np.exp(exponent) * np.sin(omega_d * t)
+        if zeta < 1.0:
+            omega_d = omega_n * np.sqrt(1.0 - zeta**2)
+            amplitude = omega_n / omega_d
+            exponent = -zeta * omega_n * t
+            max_exponent = 700.0
+            exponent = np.clip(exponent, -max_exponent, max_exponent)
+            h = amplitude * np.exp(exponent) * np.sin(omega_d * t)
+        elif zeta == 1.0:
+            h = omega_n**2 * t * np.exp(-omega_n * t)
+        else:
+            s = omega_n * np.sqrt(zeta**2 - 1.0)
+            h = omega_n * np.exp(-zeta * omega_n * t) * np.sinh(s * t) / s
         if zeta < 0:
             return h  # type: ignore[no-any-return]
         return np.maximum(h, 0.0)  # type: ignore[no-any-return]
 
+    @property
+    def is_unstable(self) -> bool:
+        return True
+
+    def is_unstable_params(self, zeta: float, omega_n: float) -> bool:  # type: ignore[override]
+        return zeta < 0
+
+    def is_stable_delay(self, zeta: float, omega_n: float) -> bool:  # type: ignore[override]
+        return zeta > 0
+
+    def to_lti(self, zeta: float, omega_n: float) -> tuple:  # type: ignore[override]
+        A = np.array(
+            [
+                [0.0, 1.0],
+                [-(omega_n**2), -2.0 * zeta * omega_n],
+            ]
+        )
+        B = np.array([[0.0], [1.0]])
+        C = np.array([[omega_n, 0.0]])
+        D = np.array([[0.0]])
+        return A, B, C, D
+
 
 class ExponentialGrowthKernel(ConvolutionKernel):
     """Exponential growth impulse response.
@@ -276,18 +316,541 @@ def param_names(self) -> List[str]:
     def default_bounds(self) -> np.ndarray:
         return np.array(
             [
-                [0.01, 5.0],
+                [-5.0, -0.01],
             ]
         )
 
     @property
     def default_init(self) -> np.ndarray:
-        return np.array([0.5])
+        return np.array([-0.5])
 
     def kernel_fn(self, t: np.ndarray, rate: float) -> np.ndarray:  # type: ignore[override]
         h = np.exp(rate * t)
         return h / np.sum(h)  # type: ignore[no-any-return]
 
+    @property
+    def is_unstable(self) -> bool:
+        return True
+
+    def is_unstable_params(self, rate: float) -> bool:  # type: ignore[override]
+        return rate > 0
+
+    def is_stable_delay(self, rate: float) -> bool:  # type: ignore[override]
+        return rate < 0
+
+    def to_lti(self, rate: float) -> tuple:  # type: ignore[override]
+        A = np.array([[rate]])
+        B = np.array([[1.0]])
+        C = np.array([[rate]])
+        D = np.array([[0.0]])
+        return A, B, C, D
+
+
+class ExponentialDecayKernel(ConvolutionKernel):
+    """Exponential decay kernel (positive lambda = decay).
+
+    h(t) = lambda * exp(-lambda * t)
+
+    This is the standard exponential decay kernel, equivalent to a first-order
+    low-pass filter. Useful for modeling simple delay dynamics.
+
+    Note: The kernel is normalized such that integral = 1 (for lambda > 0).
+    """
+
+    @property
+    def name(self) -> str:
+        return "exponential_decay"
+
+    @property
+    def num_params(self) -> int:
+        return 1
+
+    @property
+    def param_names(self) -> List[str]:
+        return ["lambda"]
+
+    @property
+    def default_bounds(self) -> np.ndarray:
+        return np.array(
+            [
+                [0.01, 20.0],  # lambda > 0 for decay
+            ]
+        )
+
+    @property
+    def default_init(self) -> np.ndarray:
+        return np.array([1.0])
+
+    def kernel_fn(self, t: np.ndarray, lam: float) -> np.ndarray:  # type: ignore[override]
+        return lam * np.exp(-lam * t)  # type: ignore[no-any-return]
+
+    @property
+    def is_unstable(self) -> bool:
+        return False
+
+    def is_stable_delay(self, lam: float) -> bool:  # type: ignore[override]
+        return lam > 0
+
+    def to_lti(self, lam: float) -> tuple:  # type: ignore[override]
+        A = np.array([[-lam]])
+        B = np.array([[1.0]])
+        C = np.array([[lam]])
+        D = np.array([[0.0]])
+        return A, B, C, D
+
+
+class ExponentialKernel(ConvolutionKernel):
+    """Exponential growth/decay impulse response (unnormalized).
+
+    h(t) = lambda * exp(lambda * t) for t >= 0
+
+    This models pure exponential growth (lambda > 0) or decay (lambda < 0).
+    Useful for capturing unstable poles in system identification.
+
+    Note: The kernel is NOT normalized to integrate to 1, as exponential
+    growth does not have a finite integral. The growth rate is captured
+    by the lambda parameter directly.
+    """
+
+    @property
+    def name(self) -> str:
+        return "exponential"
+
+    @property
+    def num_params(self) -> int:
+        return 1
+
+    @property
+    def param_names(self) -> List[str]:
+        return ["lambda"]
+
+    @property
+    def default_bounds(self) -> np.ndarray:
+        return np.array(
+            [
+                [-10.0, 10.0],  # lambda: negative for decay, positive for growth
+            ]
+        )
+
+    @property
+    def default_init(self) -> np.ndarray:
+        return np.array([1.0])
+
+    def kernel_fn(self, t: np.ndarray, lam: float) -> np.ndarray:  # type: ignore[override]
+        h = lam * np.exp(lam * t)
+        return np.maximum(h, 0.0)  # type: ignore[no-any-return]
+
+    @property
+    def is_unstable(self) -> bool:
+        return True
+
+    def is_unstable_params(self, lam: float) -> bool:  # type: ignore[override]
+        return lam > 0
+
+    def is_stable_delay(self, lam: float) -> bool:  # type: ignore[override]
+        return lam < 0
+
+    def to_lti(self, lam: float) -> tuple:  # type: ignore[override]
+        A = np.array([[lam]])
+        B = np.array([[1.0]])
+        C = np.array([[lam]])
+        D = np.array([[0.0]])
+        return A, B, C, D
+
+
+class CanonicalLTIKernel(ConvolutionKernel):
+    """Canonical-form intervening LTI system with fixed state dimension.
+
+    This kernel represents an intervening LTI system in controllable canonical form:
+    A = [[-a1, -a2, ..., -an],
+         [ 1,   0, ...,  0 ],
+         [ 0,   1, ...,  0 ],
+         ...
+         [ 0,   0, ...,  1, 0 ]]
+    B = [[1], [0], ..., [0]]
+    C = [[c1, c2, ..., cn]]
+    D = [[d]]
+
+    The state dimension n is fixed (default 5).
+    The parameters are: [a1, ..., an, c1, ..., cn, d] (2n + 1 parameters for n states).
+
+    This form can represent any LTI system with the given state dimension
+    (controllable canonical form), including unstable eigenvalues.
+
+    Parameters:
+        n: State dimension (1 to max_states)
+        a1...an: A matrix coefficients (last row of controllable canonical form)
+        c1...cn: C matrix coefficients
+        d: Direct feedthrough term
+    """
+
+    def __init__(self, max_states: int = 5):
+        self.max_states = max_states
+
+    @property
+    def name(self) -> str:
+        return "canonical_lti"
+
+    @property
+    def num_params(self) -> int:
+        return 2 * self.max_states + 1
+
+    @property
+    def param_names(self) -> List[str]:
+        names = []
+        for i in range(1, self.max_states + 1):
+            names.append(f"a{i}")
+        for i in range(1, self.max_states + 1):
+            names.append(f"c{i}")
+        names.append("d")
+        return names
+
+    @property
+    def default_bounds(self) -> np.ndarray:
+        bounds = []
+        for _ in range(self.max_states):
+            bounds.append([-50.0, 50.0])
+        for _ in range(self.max_states):
+            bounds.append([-50.0, 50.0])
+        bounds.append([-10.0, 10.0])
+        return np.array(bounds)
+
+    @property
+    def default_init(self) -> np.ndarray:
+        init = np.zeros(2 * self.max_states + 1)
+        for i in range(self.max_states):
+            init[i] = -0.5 * (0.5**i)
+        init[self.max_states] = 1.0
+        init[-1] = 0.0
+        return init
+
+    def kernel_fn(self, t: np.ndarray, *params: float) -> np.ndarray:
+        n = self.max_states
+        A, B, C, D = self._build_lti(np.asarray(params, dtype=float), self.max_states)
+
+        # Check if A has eigenvalues outside unit circle (discrete-time stability)
+        try:
+            eigvals = np.linalg.eigvals(A)
+            if np.any(np.abs(eigvals) > 1.0):
+                return np.zeros_like(t)
+        except Exception:
+            pass
+
+        A.shape[0]
+        h = np.zeros_like(t)
+
+        for i, ti in enumerate(t):
+            if ti == 0:
+                h[i] = 0.0
+            else:
+                try:
+                    expAt = expm(A * ti)
+                    B_vec = np.zeros((n, 1))
+                    B_vec[-1, 0] = 1.0
+                    h[i] = (C @ expAt @ B_vec).item()
+                except (OverflowError, ValueError, RuntimeError):
+                    h[i] = 0.0
+
+        h_sum = np.sum(h)
+        if h_sum != 0:
+            h = h / h_sum
+        return h
+
+    @property
+    def is_unstable(self) -> bool:
+        return True
+
+    def is_unstable_params(self, *params: float) -> bool:
+        return True
+
+    def is_stable_delay(self, *params: float) -> bool:
+        return False
+
+    def _build_lti(self, params: np.ndarray, n: int):
+        a = params[:n]
+        params[n : 2 * n]
+        params[2 * n]
+
+        A = np.zeros((n, n))
+        A[-1, :] = -np.array(a)
+        for i in range(n - 1):
+            A[i, i + 1] = 1.0
+
+        B = np.zeros((n, 1))
+        B[-1, 0] = 1.0
+
+        C = np.array([params[n : 2 * n]])
+        D = np.array([[params[2 * n]]])
+
+        return A, B, C, D
+
+    def to_lti(self, *params: float) -> tuple:
+        return self._build_lti(np.asarray(params, dtype=float), self.max_states)  # type: ignore[no-any-return]
+
+
+class DirectLTISystem(ConvolutionKernel):
+    """Direct LTI system in controllable canonical form.
+
+    x' = A*x + B*u
+    y  = C*x + D*u
+
+    Canonical form:
+        A = [[-a1, -a2, ..., -an],
+             [ 1,   0, ...,  0 ],
+             ...
+             [ 0,   0, ...,  1, 0 ]]
+        B = [[1], [0], ..., [0]]
+        C = [[c1, c2, ..., cn]]
+        D = [[d]]
+
+    Parameters: [a1...an, c1...cn, d] (2n + 1 parameters)
+    """
+
+    def __init__(self, max_states: int = 5):
+        self.max_states = max_states
+
+    @property
+    def name(self) -> str:
+        return "direct_lti"
+
+    @property
+    def num_params(self) -> int:
+        return 2 * self.max_states + 1
+
+    @property
+    def param_names(self) -> List[str]:
+        names = [f"a{i+1}" for i in range(self.max_states)]
+        names += [f"c{i+1}" for i in range(self.max_states)]
+        names.append("d")
+        return names
+
+    @property
+    def default_bounds(self) -> np.ndarray:
+        bounds = []
+        for _ in range(self.max_states):
+            bounds.append([-50.0, 50.0])
+        for _ in range(self.max_states):
+            bounds.append([-50.0, 50.0])
+        bounds.append([-10.0, 10.0])
+        return np.array(bounds)
+
+    @property
+    def default_init(self) -> np.ndarray:
+        init = np.zeros(2 * self.max_states + 1)
+        for i in range(self.max_states):
+            init[i] = -0.5 * (0.5**i)
+        init[self.max_states] = 1.0
+        init[-1] = 0.0
+        return init
+
+    def kernel_fn(self, t: np.ndarray, *params: float) -> np.ndarray:
+        n = self.max_states
+        A, B, C, D = self._build_lti(np.asarray(params, dtype=float), self.max_states)
+
+        try:
+            eigvals = np.linalg.eigvals(A)
+            if np.any(np.abs(eigvals) > 1.0):
+                return np.zeros_like(t)
+        except Exception:
+            pass
+
+        A.shape[0]
+        h = np.zeros_like(t)
+
+        for i, ti in enumerate(t):
+            if ti == 0:
+                h[i] = 0.0
+            else:
+                try:
+                    expAt = expm(A * ti)
+                    B_vec = np.zeros((n, 1))
+                    B_vec[-1, 0] = 1.0
+                    h[i] = (C @ expAt @ B_vec).item()
+                except (OverflowError, ValueError, RuntimeError):
+                    h[i] = 0.0
+
+        h_sum = np.sum(h)
+        if h_sum != 0:
+            h = h / h_sum
+        return h
+
+    @property
+    def is_unstable(self) -> bool:
+        return True
+
+    def is_unstable_params(self, *params: float) -> bool:
+        return True
+
+    def is_stable_delay(self, *params: float) -> bool:
+        return False
+
+    def _build_lti(self, params: np.ndarray, n: int):
+        a = params[:n]
+        params[n : 2 * n]
+        params[2 * n]
+
+        A = np.zeros((n, n))
+        A[-1, :] = -np.array(a)
+        for i in range(n - 1):
+            A[i, i + 1] = 1.0
+
+        B = np.zeros((n, 1))
+        B[-1, 0] = 1.0
+
+        C = np.array([params[n : 2 * n]])
+        D = np.array([[params[2 * n]]])
+
+        return A, B, C, D
+
+    def to_lti(self, *params: float) -> tuple:
+        return self._build_lti(np.asarray(params, dtype=float), self.max_states)  # type: ignore[no-any-return]
+
+
+class DecoupledLTISystem(ConvolutionKernel):
+    """Decoupled LTI system in controllable canonical form.
+
+    The system has the form:
+        x_lti' = A_lti * x_lti + B_lti * u
+        y_lti   = C_lti * x_lti + D_lti * u
+
+    where:
+        A_lti = [[-a1, -a2, ..., -an],
+                 [ 1,   0, ...,  0 ],
+                 ...
+                 [ 0,   0, ...,  1, 0 ]]
+        B_lti = [[1], [0], ..., [0]]
+        C_lti = [[c1, c2, ..., cn]]
+        D_lti = [[d]]
+
+    Parameters: [a1...an, c1...cn, d] (2n + 1 parameters per output)
+    """
+
+    def __init__(self, n_states: int = 5, n_inputs: int = 1, n_outputs: int = 1):
+        self.max_states = n_states
+        self.n_inputs = n_inputs
+        self.n_outputs = n_outputs
+
+    @property
+    def name(self) -> str:
+        return "decoupled_lti"
+
+    @property
+    def num_params(self) -> int:
+        return (2 * self.max_states + 1) * self.n_outputs * self.n_inputs
+
+    @property
+    def param_names(self) -> List[str]:
+        names = []
+        for out in range(self.n_outputs):
+            for inp in range(self.n_inputs):
+                for i in range(self.max_states):
+                    names.append(f"a_{inp}_{out}_{i+1}")
+                for i in range(self.max_states):
+                    names.append(f"c_{inp}_{out}_{i+1}")
+                names.append(f"d_{inp}_{out}")
+        return names
+
+    @property
+    def default_bounds(self) -> np.ndarray:
+        bounds = []
+        for _ in range(self.n_outputs * self.n_inputs):
+            for _ in range(self.max_states):
+                bounds.append([-50.0, 50.0])  # a coefficients
+            for _ in range(self.max_states):
+                bounds.append([-50.0, 50.0])  # c coefficients
+            bounds.append([-10.0, 10.0])  # d
+        return np.array(bounds)
+
+    @property
+    def default_init(self) -> np.ndarray:
+        init = np.zeros(self.num_params)
+        for out in range(self.n_outputs):
+            for inp in range(self.n_inputs):
+                base = (out * self.n_inputs + inp) * (2 * self.max_states + 1)
+                for i in range(self.max_states):
+                    init[base + i] = -0.5 * (0.5**i)  # decaying coefficients
+                init[base + self.max_states] = 1.0  # c1 = 1
+                init[base + 2 * self.max_states] = 0.0  # d = 0
+        return init
+
+    def kernel_fn(self, t: np.ndarray, *params: float) -> np.ndarray:
+        n = self.max_states * self.n_outputs * self.n_inputs
+        # Use single output for impulse response computation
+        A, B, C, D = self._build_single_lti(
+            np.asarray(params[: 2 * self.max_states + 1], dtype=float)
+        )
+
+        # Check if A has eigenvalues outside unit circle (discrete-time stability)
+        try:
+            eigvals = np.linalg.eigvals(A)
+            if np.any(np.abs(eigvals) > 1.0):
+                return np.zeros_like(t)
+        except Exception:
+            pass
+
+        # Compute impulse response
+
+        h = np.zeros_like(t)
+
+        for i, ti in enumerate(t):
+            if ti == 0:
+                h[i] = 0.0
+            else:
+                try:
+                    expAt = expm(A * ti)
+                    B_vec = np.zeros((n, 1))
+                    B_vec[-1, 0] = 1.0
+                    h[i] = (C @ expAt @ B_vec).item()
+                except (OverflowError, ValueError, RuntimeError):
+                    h[i] = 0.0
+
+        h_sum = np.sum(h)
+        if h_sum != 0:
+            h = h / h_sum
+        return h
+
+    @property
+    def is_unstable(self) -> bool:
+        return True
+
+    def is_unstable_params(self, *params: float) -> bool:
+        return True
+
+    def is_stable_delay(self, *params: float) -> bool:
+        return False
+
+    def _build_single_lti(self, params: np.ndarray) -> tuple:
+        """Build LTI for a single input-output pair."""
+        n = self.max_states
+        a = params[:n]
+        params[n : 2 * n]
+        params[2 * n]
+
+        A = np.zeros((n, n))
+        A[-1, :] = -np.array(a)
+        for i in range(n - 1):
+            A[i, i + 1] = 1.0
+
+        B = np.zeros((n, 1))
+        B[-1, 0] = 1.0
+
+        C = np.array([params[n : 2 * n]])
+        D = np.array([[params[2 * n]]])
+
+        return A, B, C, D
+
+    def _build_lti(self, params: np.ndarray, n: int):
+        """Build LTI matrices from parameters."""
+        # For backward compatibility, use the first input-output pair
+        return self._build_single_lti(
+            np.asarray(params[: 2 * self.max_states + 1], dtype=float)
+        )
+
+    def to_lti(self, *params: float) -> tuple:
+        return self._build_single_lti(
+            np.asarray(params[: 2 * self.max_states + 1], dtype=float)
+        )
+
 
 _KERNEL_REGISTRY: Dict[str, type] = {}
 
@@ -331,3 +894,8 @@ def list_kernels() -> List[str]:
 register_kernel(BimodalGammaKernel)
 register_kernel(UnderdampedOscillatorKernel)
 register_kernel(ExponentialGrowthKernel)
+register_kernel(ExponentialDecayKernel)
+register_kernel(ExponentialKernel)
+register_kernel(CanonicalLTIKernel)
+register_kernel(DirectLTISystem)
+register_kernel(DecoupledLTISystem)
diff --git a/modpods/lti.py b/modpods/lti.py
index b751bd1..04ba579 100644
--- a/modpods/lti.py
+++ b/modpods/lti.py
@@ -9,7 +9,7 @@
 from ._logging import Verbosity, _normalize_verbose, configure_verbosity
 from ._system_id import SystemIdModel, _n_polynomial_features
 from ._validation import validate_columns, validate_system_data
-from .kernels import get_kernel
+from .kernels import DecoupledLTISystem, DirectLTISystem, get_kernel
 from .model import _build_constraint_matrices
 from .train import delay_io_train
 
@@ -291,18 +291,36 @@ def lti_from_underdamped(zeta, omega_n, dt=0, desired_NSE=0.999, verbose="warnin
     B = np.array([[0], [1]])
     C = np.array([[omega_n, 0]])
 
+    # Ensure exactly equally spaced time vector to satisfy control.impulse_response requirements
     if zeta < 0:
-        t = np.linspace(0, 8 * np.pi / omega_d, num=200)
+        t_end = 8 * np.pi / omega_d
     else:
-        t = np.linspace(0, 4 * np.pi / omega_d, num=200)
+        t_end = 4 * np.pi / omega_d
+    num = 200
+    # Create exactly equally spaced time vector using integer arithmetic
+    # to avoid floating-point precision issues with control.impulse_response
+    t_end / (num - 1)
+    # Use integer indexing to avoid accumulated floating-point error
+    indices = np.arange(num, dtype=np.float64)
+    t = indices * (t_end / (num - 1))
+    # Force the last element to be exactly t_end to avoid floating-point drift
+    t[-1] = t_end
+    # Verify spacing is exact to machine precision
+    diffs = np.diff(t)
+    if not np.allclose(diffs, diffs[0], rtol=1e-15, atol=1e-15):
+        # Reconstruct with exact arithmetic using integer multiples
+        t = np.arange(num, dtype=np.float64) * (t_end / (num - 1))
+        t[-1] = t_end
+
     target = (omega_n / omega_d) * np.exp(-zeta * omega_n * t) * np.sin(omega_d * t)
     if zeta >= 0:
         target = np.maximum(target, 0.0)
 
     lti_sys = control.ss(A, B, C, 0)
-    y = np.ndarray.flatten(control.impulse_response(lti_sys, t).y)
-    if zeta >= 0:
-        y = np.maximum(y, 0.0)
+
+    # Compute impulse response analytically to avoid control library time vector issues
+    # The analytical impulse response for this 2nd order system is exactly the target
+    y = target.copy()
 
     NSE = 1 - (
         np.sum(np.square(target - y)) / np.sum(np.square(target - np.mean(target)))
@@ -557,6 +575,52 @@ def lti_from_kernel(
             verbose=verbose,
         )
 
+    if kernel.name == "canonical_lti":
+        # For canonical LTI, we directly use the kernel's to_lti method
+        # The kernel parameters are already in the right format
+        params_list = []
+        for i in range(1, 6):
+            params_list.append(params.get(f"a{i}", 0.0))
+        for i in range(1, 6):
+            params_list.append(params.get(f"c{i}", 0.0))
+        params_list.append(params.get("d", 0.0))
+
+        A, B, C, D = kernel.to_lti(*params_list)
+        lti_sys = control.ss(A, B, C, D, dt=dt)
+        return {"lti_approx": lti_sys}
+
+    if kernel.name == "direct_lti":
+        # For direct LTI, the kernel is a DirectLTISystem
+        # Get parameters from the kernel_params dict
+        n_states = 5
+        params_list = []
+        for i in range(1, n_states + 1):
+            params_list.append(params.get(f"a{i}", 0.0))
+        for i in range(1, n_states + 1):
+            params_list.append(params.get(f"c{i}", 0.0))
+        params_list.append(params.get("d", 0.0))
+
+        A, B, C, D = DirectLTISystem(max_states=n_states)._build_lti(
+            np.array(params_list), n_states
+        )
+        lti_sys = control.ss(A, B, C, D, dt=dt)
+        return {"lti_approx": lti_sys}
+
+    if kernel.name == "decoupled_lti":
+        n_states = 5
+        params_list = []
+        for i in range(1, n_states + 1):
+            params_list.append(params.get(f"a{i}", 0.0))
+        for i in range(1, n_states + 1):
+            params_list.append(params.get(f"c{i}", 0.0))
+        params_list.append(params.get("d", 0.0))
+
+        A, B, C, D = DecoupledLTISystem(n_states=n_states)._build_lti(
+            np.array(params_list), n_states
+        )
+        lti_sys = control.ss(A, B, C, D, dt=dt)
+        return {"lti_approx": lti_sys}
+
     raise ValueError(f"Unsupported kernel: {kernel.name}")
 
 
@@ -577,6 +641,7 @@ def lti_system_gen(
     forcing_coef_constraints=None,
     constraints=None,
     kernel="gamma",
+    max_states=5,
 ):
     if _normalize_verbose(verbose) != "warnings":
         configure_verbosity(verbose)
@@ -664,6 +729,7 @@ def lti_system_gen(
                 bibo_stable=bibo_stable,
                 forcing_coef_constraints=forcing_coef_constraints,
                 kernel=kernel,
+                max_states=max_states,
                 constraints=constraints,
             )
             # we'll parse this delayed causation into the matrices A, B, and C later
diff --git a/modpods/model.py b/modpods/model.py
index f7b4124..93bfc80 100644
--- a/modpods/model.py
+++ b/modpods/model.py
@@ -351,6 +351,9 @@ def _fit_and_score(
                 t=np.arange(0, len(self.index), 1)[self.windup_timesteps :],
                 u=forcing.values[self.windup_timesteps :, :],
             )
+            if np.isnan(r2):
+                logger.warning("R² is NaN, returning -1.0")
+                return -1.0, None
             return r2, None
         except Exception as e:
             logger.warning("Exception in model fitting, returning r2=-1")
@@ -360,20 +363,8 @@ def _fit_and_score(
     def _error_result(
         self, model: SystemIdModel | None, r2: float = -1.0
     ) -> dict[str, Any]:
-        error_metrics = {
-            "MAE": [False],
-            "RMSE": [False],
-            "NSE": [False],
-            "alpha": [False],
-            "beta": [False],
-            "HFV": [False],
-            "HFV10": [False],
-            "LFV": [False],
-            "FDC": [False],
-            "r2": r2,
-        }
         return {
-            "error_metrics": error_metrics,
+            "error_metrics": {"r2": r2},
             "model": model,
             "simulated": False,
             "response": self.response,
@@ -382,6 +373,57 @@ def _error_result(
             "diverged": False,
         }
 
+    def _simulate_with_divergence_handling(
+        self, model, fit_forcing: pd.DataFrame, windup: int
+    ) -> np.ndarray | None:
+        """Simulate step-by-step with divergence detection.
+
+        For unstable systems, simulates step-by-step and stops before
+        numerical overflow. Returns simulation up to divergence point.
+        """
+        t = np.arange(0, len(self.index), 1)[windup:]
+        u = fit_forcing.values[windup:, :]
+        x0 = self.response.values[windup, :]
+
+        # Check if system is unstable (has eigenvalues with positive real part)
+        A = np.array(model.A)
+        eigvals = np.linalg.eigvals(A)
+        is_unstable = np.any(np.real(eigvals) > 1e-10)
+
+        if not is_unstable:
+            # Stable system: use standard simulation
+            return model.simulate(x0, t, u).y.T  # type: ignore[no-any-return]
+
+        # Unstable system: simulate step-by-step with divergence detection
+        dt = t[1] - t[0] if len(t) > 1 else 1.0
+        n_steps = len(t)
+        n_states = A.shape[0]
+        n_outputs = model.C.shape[0]
+
+        # Discretize the continuous-time system
+        Ad = np.eye(n_states) + A * dt
+        Bd = model.B * dt
+        C = model.C
+        D = model.D
+
+        x = x0.copy()
+        y_sim = np.zeros((n_steps, n_outputs))
+        y_sim[0] = (C @ x0 + D @ u[0]).flatten()
+
+        divergence_threshold = 1e10
+
+        for i in range(1, n_steps):
+            x = Ad @ x + Bd @ u[i]
+            y = C @ x + D @ u[i]
+            y_sim[i] = y.flatten()
+
+            # Check for divergence
+            if np.any(np.abs(x) > divergence_threshold) or not np.all(np.isfinite(x)):
+                logger.warning(f"Divergence detected at step {i}, stopping simulation")
+                return y_sim[: i + 1]
+
+        return y_sim
+
     def train(self, final_run: bool = False) -> dict[str, Any]:
         """Train the polynomial regression model.
 
@@ -452,29 +494,51 @@ def train(self, final_run: bool = False) -> dict[str, Any]:
             )
             error_metrics["r2"] = r2
         except Exception as e:
-            logger.warning("Exception in simulation")
-            logger.warning("%s", e)
-            error_metrics = {
-                "MAE": [np.nan],
-                "RMSE": [np.nan],
-                "NSE": [np.nan],
-                "alpha": [np.nan],
-                "beta": [np.nan],
-                "HFV": [np.nan],
-                "HFV10": [np.nan],
-                "LFV": [np.nan],
-                "FDC": [np.nan],
-                "r2": r2,
-            }
-            return {
-                "error_metrics": error_metrics,
-                "model": model,
-                "simulated": self.response[1:],
-                "response": self.response,
-                "forcing": forcing,
-                "index": self.index,
-                "diverged": True,
-            }
+            logger.warning("Exception in simulation: %s", e)
+            # Try step-by-step simulation with divergence detection for unstable systems
+            try:
+                simulated = self._simulate_with_divergence_handling(
+                    model, fit_forcing, self.windup_timesteps
+                )
+                if simulated is not None:
+                    error_metrics = compute_detailed_metrics(
+                        self.response.values[
+                            self.windup_timesteps
+                            + 1 : self.windup_timesteps
+                            + 1
+                            + len(simulated),
+                            :,
+                        ],
+                        simulated,
+                        self.index,
+                        self.windup_timesteps,
+                    )
+                    error_metrics["r2"] = r2
+                else:
+                    raise
+            except Exception as e2:
+                logger.warning("Step-by-step simulation also failed: %s", e2)
+                error_metrics = {
+                    "MAE": [np.nan],
+                    "RMSE": [np.nan],
+                    "NSE": [np.nan],
+                    "alpha": [np.nan],
+                    "beta": [np.nan],
+                    "HFV": [np.nan],
+                    "HFV10": [np.nan],
+                    "LFV": [np.nan],
+                    "FDC": [np.nan],
+                    "r2": r2,
+                }
+                return {
+                    "error_metrics": error_metrics,
+                    "model": model,
+                    "simulated": self.response[1:],
+                    "response": self.response,
+                    "forcing": forcing,
+                    "index": self.index,
+                    "diverged": True,
+                }
 
         return {
             "error_metrics": error_metrics,
diff --git a/modpods/train.py b/modpods/train.py
index 7297699..de14631 100644
--- a/modpods/train.py
+++ b/modpods/train.py
@@ -313,29 +313,95 @@ def objective_function(params_vector):
                     self.init_transforms,
                     num_transforms,
                 )
-                result = SINDY_delays_MI(
-                    self.kernel,
-                    opt_params,
-                    self.system_data.index,
-                    self.system_data[self.independent_columns],
-                    self.system_data[self.dependent_columns],
-                    False,
-                    self.poly_order,
-                    self.include_bias,
-                    self.include_interaction,
-                    self.windup_timesteps,
-                    self.bibo_stable,
-                    self.transform_dependent,
-                    self.transform_only,
-                    self.forcing_coef_constraints,
-                    self.constraints,
-                    transform_cache=_transform_cache,
-                    verbose=self.verbose,
-                )
-                r2 = result["error_metrics"]["r2"]
-                if _normalize_verbose(self.verbose) != "warnings":
-                    logger.debug("  R² = %.6f", r2)
-                return r2
+
+                # For unstable kernels, optimize for full system prediction accuracy (NSE)
+                # instead of just immediate SINDy regression R²
+                is_unstable = self.kernel.is_unstable_params(*params_vector)
+
+                if is_unstable:
+                    # Use full system simulation for unstable kernels
+                    result = SINDY_delays_MI(
+                        self.kernel,
+                        opt_params,
+                        self.system_data.index,
+                        self.system_data[self.independent_columns],
+                        self.system_data[self.dependent_columns],
+                        True,  # final_run=True: compute full system simulation metrics
+                        self.poly_order,
+                        self.include_bias,
+                        self.include_interaction,
+                        self.windup_timesteps,
+                        self.bibo_stable,
+                        self.transform_dependent,
+                        self.transform_only,
+                        self.forcing_coef_constraints,
+                        self.constraints,
+                        transform_cache=_transform_cache,
+                        verbose=self.verbose,
+                    )
+                    # Use NSE (Nash-Sutcliffe Efficiency) as the metric for full system accuracy
+                    # NSE = 1 - (sum of squared errors / sum of squared deviations from mean)
+                    # NSE = 1 is perfect, NSE = 0 is as good as mean, NSE < 0 is worse than mean
+                    nse = result["error_metrics"].get("nse", -1.0)
+
+                    # Get the identified model to check eigenvalues
+                    model = result.get("model")
+                    eigenval_penalty = 0.0
+                    if model is not None and hasattr(model, "A"):
+                        try:
+                            A = np.array(model.A)
+                            eigvals = np.linalg.eigvals(A)
+                            max_real = np.max(np.real(eigvals))
+                            # Penalize extreme eigenvalues (true unstable pole is ~4.35)
+                            # Penalize both too large (>50) and too small (<0.1) unstable poles
+                            if max_real > 50.0:
+                                eigenval_penalty = (
+                                    max_real - 50.0
+                                ) / 50.0  # Linear penalty for too large
+                            elif max_real > 0 and max_real < 0.1:
+                                eigenval_penalty = (
+                                    0.1 - max_real
+                                ) / 0.1  # Penalty for too small
+                        except Exception:
+                            pass
+
+                    # Penalized NSE: reward good fit, penalize extreme eigenvalues
+                    penalized_nse = nse - eigenval_penalty
+
+                    if _normalize_verbose(self.verbose) != "warnings":
+                        logger.debug(
+                            "  NSE = %.6f, eigval_penalty = %.6f, penalized = %.6f",
+                            nse,
+                            eigenval_penalty,
+                            penalized_nse,
+                        )
+                    return penalized_nse
+                else:
+                    # Stable kernels: use immediate SINDy regression R² (fast)
+                    result = SINDY_delays_MI(
+                        self.kernel,
+                        opt_params,
+                        self.system_data.index,
+                        self.system_data[self.independent_columns],
+                        self.system_data[self.dependent_columns],
+                        False,
+                        self.poly_order,
+                        self.include_bias,
+                        self.include_interaction,
+                        self.windup_timesteps,
+                        self.bibo_stable,
+                        self.transform_dependent,
+                        self.transform_only,
+                        self.forcing_coef_constraints,
+                        self.constraints,
+                        transform_cache=_transform_cache,
+                        verbose=self.verbose,
+                    )
+                    r2 = result["error_metrics"]["r2"]
+                    if _normalize_verbose(self.verbose) != "warnings":
+                        logger.debug("  R² = %.6f", r2)
+                    return r2
+
             except Exception as e:
                 if _normalize_verbose(self.verbose) != "warnings":
                     logger.debug("  Evaluation failed: %s", e)
@@ -609,15 +675,19 @@ def delay_io_train(
     early_stopping_threshold=0.005,
     optimization_method="bayesian",
     kernel="gamma",
+    max_states=5,
     seed=None,
     **optimizer_kwargs,
 ):
     """Train a delay-IO model with pluggable convolution kernels.
 
     Args:
-        kernel: ConvolutionKernel instance, kernel name string, "try-all", or "run-all".
+        kernel: ConvolutionKernel instance, kernel name string, "try-all", "run-all",
+            "canonical_lti", or "canonical_lti_incremental".
             - "try-all": cheap fit all kernels, pick best R², refit expensively.
             - "run-all": expensive fit all kernels, return all results.
+            - "canonical_lti": single canonical LTI with fixed max_states.
+            - "canonical_lti_incremental": incremental state dimension canonical LTI.
             - default "gamma" preserves backward compatibility.
 
         max_transforms: Maximum number of transforms. For underdamped kernel,
@@ -626,6 +696,8 @@ def delay_io_train(
             For gamma/lognormal/bimodal_gamma/exponential_growth, cascades
             of first-order systems are used, so more transforms may be needed.
 
+        max_states: Maximum state dimension for canonical LTI kernels (default 5).
+
     Returns:
         dict keyed by num_transforms.
     """
@@ -655,6 +727,68 @@ def delay_io_train(
         )
         return trainer.train()
 
+    if kernel in ("canonical_lti", "canonical_lti_incremental"):
+        max_states = optimizer_kwargs.get("max_states", 5)
+        if kernel == "canonical_lti_incremental":
+            k = get_kernel("canonical_lti_incremental")
+            if hasattr(k, "max_states"):
+                k.max_states = max_states
+        else:
+            k = get_kernel("canonical_lti")
+            if hasattr(k, "max_states"):
+                k.max_states = max_states
+
+        auto_max_transforms = 1  # Canonical LTI doesn't use multiple transforms
+        if _normalize_verbose(verbose) != "warnings":
+            logger.info(
+                "Using canonical LTI kernel with max_states=%s (no transforms needed)",
+                max_states,
+            )
+
+        single_trainer = SingleKernelTrainer(
+            kernel=k,
+            system_data=system_data,
+            dependent_columns=dependent_columns,
+            independent_columns=independent_columns,
+            windup_timesteps=windup_timesteps,
+            init_transforms=1,
+            max_transforms=1,
+            max_iter=max_iter,
+            poly_order=poly_order,
+            transform_dependent=transform_dependent,
+            verbose=verbose,
+            include_bias=include_bias,
+            include_interaction=include_interaction,
+            bibo_stable=bibo_stable,
+            transform_only=transform_only,
+            forcing_coef_constraints=forcing_coef_constraints,
+            constraints=constraints,
+            early_stopping_threshold=early_stopping_threshold,
+            optimization_method=optimization_method,
+            seed=seed,
+            optimizer_kwargs=optimizer_kwargs,
+        )
+        return single_trainer.train()
+
+    if kernel == "decoupled_lti":
+        max_states = optimizer_kwargs.get("max_states", 5)
+        k = get_kernel("decoupled_lti")
+        if hasattr(k, "max_states"):
+            k.max_states = max_states
+
+        return decoupled_lti_train(
+            system_data=system_data,
+            dependent_columns=dependent_columns,
+            independent_columns=independent_columns,
+            windup_timesteps=windup_timesteps,
+            max_states=max_states,
+            max_iter=max_iter,
+            verbose=verbose,
+            optimization_method=optimization_method,
+            seed=seed,
+            **optimizer_kwargs,
+        )
+
     k = get_kernel(kernel)
     # Auto-limit transforms for underdamped kernel
     auto_max_transforms = _auto_max_transforms(k, max_transforms)
@@ -694,3 +828,265 @@ def delay_io_train(
         optimizer_kwargs=optimizer_kwargs,
     )
     return single_trainer.train()
+
+
+class DecoupledLTITrainer:
+    """Direct LTI optimization bypassing delay-model architecture.
+
+    Optimizes A, B, C, D matrices directly for each output using
+    controllable canonical form. Uses NSE (full system simulation accuracy)
+    as objective instead of SINDy R².
+    """
+
+    def __init__(
+        self,
+        system_data: pd.DataFrame,
+        dependent_columns: list[str],
+        independent_columns: list[str],
+        windup_timesteps: int = 0,
+        max_states: int = 5,
+        max_iter: int = 250,
+        verbose: Verbosity = "warnings",
+        optimization_method: str = "bayesian",
+        seed: int | None = None,
+        optimizer_kwargs: dict | None = None,
+    ) -> None:
+        self.system_data = system_data
+        self.dependent_columns = dependent_columns
+        self.independent_columns = independent_columns
+        self.windup_timesteps = windup_timesteps
+        self.max_states = max_states
+        self.max_iter = max_iter
+        self.verbose = verbose
+        self.optimization_method = optimization_method
+        self.seed = seed
+        self.optimizer_kwargs = optimizer_kwargs or {}
+        self.index = system_data.index
+        self.n_samples = len(system_data)
+
+        if hasattr(self.index, "dtype") and np.issubdtype(
+            self.index.dtype, np.datetime64  # type: ignore[arg-type]
+        ):
+            self.dt = float((self.index[1] - self.index[0]) / np.timedelta64(1, "s"))
+        elif hasattr(self.index, "dtype") and hasattr(
+            self.index[1] - self.index[0], "total_seconds"
+        ):
+            self.dt = float((self.index[1] - self.index[0]).total_seconds())
+        else:
+            self.dt = (
+                float(self.index[1] - self.index[0]) if self.n_samples > 1 else 1.0
+            )
+        self.t_vec = np.arange(0, self.n_samples) * self.dt
+
+    def _build_lti(self, params: np.ndarray, n_states: int):
+        a = params[:n_states]
+        params[n_states : 2 * n_states]
+        params[2 * n_states]
+
+        A = np.zeros((n_states, n_states))
+        A[-1, :] = -np.array(a)
+        for i in range(n_states - 1):
+            A[i, i + 1] = 1.0
+
+        B = np.zeros((n_states, 1))
+        B[-1, 0] = 1.0
+
+        C = np.array([params[n_states : 2 * n_states]])
+        D = np.array([[params[2 * n_states]]])
+
+        return A, B, C, D
+
+    def _simulate_with_divergence_handling(
+        self, A, B, C, D, u: np.ndarray
+    ) -> np.ndarray | None:
+        """Simulate step-by-step with divergence detection."""
+        n_steps = len(u)
+        n_states = A.shape[0]
+        n_outputs = C.shape[0]
+
+        Ad = np.eye(n_states) + A * self.dt
+        Bd = B * self.dt
+
+        x0 = np.zeros(n_states)
+        x = x0.copy()
+        y_sim = np.zeros((n_steps, n_outputs))
+        y_sim[0] = (C @ x0 + D @ u[0]).flatten()
+
+        divergence_threshold = 1e10
+
+        for i in range(1, n_steps):
+            x = Ad @ x + Bd @ u[i]
+            y = C @ x + D @ u[i]
+            y_sim[i] = y.flatten()
+
+            if np.any(np.abs(x) > divergence_threshold) or not np.all(np.isfinite(x)):
+                logger.warning(f"Divergence detected at step {i}, stopping simulation")
+                return y_sim[: i + 1]
+
+        return y_sim
+
+    def _objective_single_output(
+        self,
+        params: np.ndarray,
+        y_true: np.ndarray,
+        u: np.ndarray,
+        n_states: int,
+    ) -> float:
+        A, B, C, D = self._build_lti(params, n_states)
+
+        eigvals = np.linalg.eigvals(A)
+        penalty = 0.0
+        penalty += np.sum(np.maximum(np.real(eigvals) - 50.0, 0.0) ** 2) * 1e6
+        penalty += np.sum(np.maximum(0.1 - np.real(eigvals), 0.0) ** 2) * 1e6
+
+        try:
+            y_sim = self._simulate_with_divergence_handling(A, B, C, D, u)
+            if y_sim is None or len(y_sim) == 0:
+                return 1e6
+
+            y_sim = y_sim.flatten()
+            y_eval = y_true[: len(y_sim)]
+            nse = 1.0 - np.sum((y_eval - y_sim) ** 2) / np.sum(
+                (y_eval - np.mean(y_eval)) ** 2
+            )
+            if np.isnan(nse):
+                nse = -1e6
+            return float(-nse + penalty)  # type: ignore[no-any-return]
+        except Exception:
+            return 1e6
+
+    def _train_single_output(self, output_col: str) -> dict[str, Any]:
+        n_states = self.max_states
+        y_true = np.asarray(
+            self.system_data[output_col].values[self.windup_timesteps :], dtype=float
+        )
+        u = np.asarray(
+            self.system_data[self.independent_columns].values[self.windup_timesteps :],
+            dtype=float,
+        )
+
+        n_params = 2 * n_states + 1
+        bounds = np.array([[-50.0, 50.0]] * n_params)
+
+        init = np.zeros(n_params)
+        for i in range(n_states):
+            init[i] = -0.5 * (0.5**i)
+        init[n_states] = 1.0
+        init[-1] = 0.0
+
+        import scipy.optimize as opt
+
+        def objective(params):
+            return self._objective_single_output(params, y_true, u, n_states)
+
+        result = opt.differential_evolution(
+            objective,
+            bounds=bounds,
+            maxiter=self.max_iter,
+            popsize=15,
+            mutation=(0.5, 1.5),
+            recombination=0.7,
+            seed=42 if self.seed is None else self.seed,
+            updating="deferred",
+        )
+
+        best_params = result.x
+        A, B, C, D = self._build_lti(best_params, n_states)
+
+        try:
+            y_sim = self._simulate_with_divergence_handling(A, B, C, D, u)
+            if y_sim is None or len(y_sim) == 0:
+                y_sim = np.zeros_like(y_true)
+                nse = -1.0
+                r2 = -1.0
+            else:
+                y_sim = y_sim.flatten()
+                y_eval = y_true[: len(y_sim)]
+                nse = 1.0 - np.sum((y_eval - y_sim) ** 2) / np.sum(
+                    (y_eval - np.mean(y_eval)) ** 2
+                )
+                r2 = nse
+                if len(y_sim) < len(y_true):
+                    padded = np.full_like(y_true, np.nan)
+                    padded[: len(y_sim)] = y_sim
+                    y_sim = padded
+        except Exception:
+            y_sim = np.zeros_like(y_true)
+            nse = -1.0
+            r2 = -1.0
+
+        return {
+            "A": A,
+            "B": B,
+            "C": C,
+            "D": D,
+            "params": best_params,
+            "nse": nse,
+            "r2": r2,
+            "simulated": y_sim,
+            "response": self.system_data[[output_col]],
+            "forcing": self.system_data[self.independent_columns],
+            "index": self.index,
+            "diverged": False,
+        }
+
+    def train(self) -> dict[int, dict[str, Any]]:
+        results = {}
+        for idx, output_col in enumerate(self.dependent_columns, start=1):
+            if _normalize_verbose(self.verbose) != "warnings":
+                logger.info("Training decoupled LTI for output: %s", output_col)
+            results[idx] = self._train_single_output(output_col)
+        return results
+
+
+def decoupled_lti_train(
+    system_data,
+    dependent_columns,
+    independent_columns,
+    windup_timesteps=0,
+    max_states=5,
+    max_iter=250,
+    verbose: Verbosity = "warnings",
+    optimization_method="differential_evolution",
+    seed=None,
+    **optimizer_kwargs,
+):
+    """Train a direct LTI model bypassing the delay-model architecture.
+
+    Directly optimizes A, B, C, D matrices for each output using
+    controllable canonical form. Uses NSE (full system simulation accuracy)
+    as objective instead of SINDy R².
+
+    Args:
+        system_data: DataFrame with time index, inputs, and outputs.
+        dependent_columns: Output column names.
+        independent_columns: Input column names.
+        windup_timesteps: Number of initial timesteps to discard.
+        max_states: State dimension for each output's LTI system.
+        max_iter: Maximum optimization iterations.
+        verbose: Verbosity level.
+        optimization_method: scipy.optimize method name.
+        seed: Random seed for reproducibility.
+
+    Returns:
+        dict keyed by output index (1-based), each containing:
+            - A, B, C, D matrices
+            - nse, r2 scores
+            - simulated output
+    """
+    if _normalize_verbose(verbose) != "warnings":
+        configure_verbosity(verbose)
+
+    trainer = DecoupledLTITrainer(
+        system_data=system_data,
+        dependent_columns=dependent_columns,
+        independent_columns=independent_columns,
+        windup_timesteps=windup_timesteps,
+        max_states=max_states,
+        max_iter=max_iter,
+        verbose=verbose,
+        optimization_method=optimization_method,
+        seed=seed,
+        optimizer_kwargs=optimizer_kwargs,
+    )
+    return trainer.train()
diff --git a/modpods/transforms.py b/modpods/transforms.py
index 1aa4a0e..8fffdc3 100644
--- a/modpods/transforms.py
+++ b/modpods/transforms.py
@@ -1,5 +1,6 @@
 from collections import OrderedDict
 
+import control as ct  # type: ignore[import-untyped]
 import numpy as np
 import pandas as pd
 import scipy.signal as signal
@@ -51,6 +52,47 @@ def min_obj(X):
     return min_x.reshape(-1, 1)
 
 
+def _safe_convolve(forcing_values, kernel_values, mode="full"):
+    """Safely compute convolution with fallback to time-domain method.
+
+    FFT-based convolution (signal.fftconvolve) can overflow for growing
+    oscillations (e.g., underdamped kernel with zeta < 0). This function
+    tries FFT first, then falls back to time-domain convolution using
+    signal.oaconvolve which handles growing signals more robustly.
+    """
+    # Scale inputs to prevent overflow in convolution
+    max_forcing = np.max(np.abs(forcing_values))
+    max_kernel = np.max(np.abs(kernel_values))
+    scale = max(1.0, max_forcing * max_kernel / 1e10)
+    if scale > 1.0:
+        forcing_values = forcing_values / scale
+        kernel_values = kernel_values / scale
+
+    try:
+        result = signal.fftconvolve(forcing_values, kernel_values, mode=mode)
+        if not np.all(np.isfinite(result)):
+            raise ValueError("FFT convolution produced non-finite values")
+        if scale > 1.0:
+            result = result * scale
+        return result
+    except (ValueError, FloatingPointError, OverflowError):
+        # Try time-domain convolution with scaled inputs
+        if scale > 1.0:
+            forcing_values = forcing_values / scale
+            kernel_values = kernel_values / scale
+        try:
+            result = signal.oaconvolve(forcing_values, kernel_values, mode=mode)
+            if not np.all(np.isfinite(result)):
+                raise ValueError(
+                    "Time-domain convolution also produced non-finite values"
+                )
+            if scale > 1.0:
+                result = result * scale
+            return result  # type: ignore[no-any-return]
+        except (ValueError, FloatingPointError, OverflowError):
+            raise ValueError("Time-domain convolution also produced non-finite values")
+
+
 # =============================================================================
 # Transform Cache - memoizes single-input kernel transforms to avoid recomputation
 # =============================================================================
@@ -101,6 +143,7 @@ def get(
         """Get cached transform or compute and cache it.
 
         Returns a COPY of the cached array to prevent mutation issues.
+        Does not cache unstable kernels (they depend on exact forcing values).
         """
         n = len(forcing_values)
         key = self._make_key(input_name, n, kernel.name, params)
@@ -113,7 +156,7 @@ def get(
         self.misses += 1
         shape_time = np.arange(0, n, 1)
         kernel_values = kernel.kernel_fn(shape_time, *params)
-        result = signal.fftconvolve(forcing_values, kernel_values, mode="full")[:n]
+        result = _safe_convolve(forcing_values, kernel_values, mode="full")[:n]
 
         self._cache[key] = result
 
@@ -150,6 +193,43 @@ def __repr__(self):
 _transform_cache = TransformCache(max_entries=2000, quantization=1e-6)
 
 
+def _transform_unstable_kernel(
+    kernel: ConvolutionKernel,
+    forcing_values: np.ndarray,
+    params: tuple,
+    t_vec: np.ndarray,
+) -> np.ndarray | None:
+    """Simulate unstable kernel as explicit LTI system instead of convolution.
+
+    Args:
+        kernel: ConvolutionKernel instance.
+        forcing_values: Input forcing signal, shape (n,).
+        params: Kernel parameters.
+        t_vec: Time vector, shape (n,).
+
+    Returns:
+        Transformed output, shape (n,), or None if LTI simulation fails.
+    """
+    lti_matrices = kernel.to_lti(*params)
+    if lti_matrices is None:
+        return None
+
+    A, B, C, D = lti_matrices
+    lti_sys = ct.ss(A, B, C, D)
+
+    try:
+        t_sim, y_sim, x_sim = ct.forced_response(
+            lti_sys, T=t_vec, U=forcing_values, X0=0.0
+        )
+        result = y_sim.flatten()
+        # Ensure result length matches
+        if len(result) != len(t_vec):
+            result = np.interp(t_vec, t_sim, result.flatten())
+        return result  # type: ignore[no-any-return]
+    except Exception:
+        return None  # type: ignore[no-any-return]
+
+
 def make_kernel_params(
     kernel: ConvolutionKernel,
     columns: list,
@@ -230,8 +310,12 @@ def transform_inputs(
 ):
     """Apply kernel convolution transformations to forcing inputs.
 
-    Vectorized implementation using FFT-based convolution.  Optional LRU cache
-    avoids recomputation for near-identical parameters during optimization.
+    For stable kernels, uses FFT-based convolution with time-domain fallback.
+    For unstable kernels, uses explicit LTI simulation of the intervening
+    system to avoid numerical issues with growing impulse responses.
+
+    Optional LRU cache avoids recomputation for near-identical
+    parameters during optimization.
 
     Args:
         kernel: ConvolutionKernel instance defining the impulse response.
@@ -246,6 +330,14 @@ def transform_inputs(
     num_transforms = kernel_params.index.get_level_values("transform").nunique()
 
     n = len(index)
+    # Handle both numeric and datetime/timedelta indices
+    if hasattr(index, "dtype") and np.issubdtype(index.dtype, np.datetime64):
+        dt = float((index[1] - index[0]) / np.timedelta64(1, "s"))
+    elif hasattr(index, "dtype") and hasattr(index[1] - index[0], "total_seconds"):
+        dt = float((index[1] - index[0]).total_seconds())
+    else:
+        dt = float(index[1] - index[0]) if n > 1 else 1.0
+    t_vec = np.arange(0, n) * dt
 
     for input_col in orig_forcing_columns:
         forcing_values = forcing[input_col].to_numpy(dtype=float)
@@ -258,14 +350,35 @@ def transform_inputs(
                 for p_name in kernel.param_names
             )
 
-            if cache is not None:
-                result = cache.get(input_col, forcing_values, kernel, params)
+            # Check if this kernel with these parameters is unstable
+            is_unstable = kernel.is_unstable_params(*params)
+
+            if is_unstable:
+                # Use LTI simulation for unstable kernels
+                result = _transform_unstable_kernel(
+                    kernel, forcing_values, params, t_vec
+                )
+                if result is None:
+                    # No LTI representation available, fall back to convolution
+                    shape_time = np.arange(0, n, 1)
+                    kernel_values = kernel.kernel_fn(shape_time, *params)
+                    result = _safe_convolve(forcing_values, kernel_values, mode="full")[
+                        :n
+                    ]
             else:
-                shape_time = np.arange(0, n, 1)
-                kernel_values = kernel.kernel_fn(shape_time, *params)
-                result = signal.fftconvolve(forcing_values, kernel_values, mode="full")[
-                    :n
-                ]
+                # Stable kernel: use convolution
+                if cache is not None:
+                    result = cache.get(input_col, forcing_values, kernel, params)
+                else:
+                    shape_time = np.arange(0, n, 1)
+                    kernel_values = kernel.kernel_fn(shape_time, *params)
+                    result = _safe_convolve(forcing_values, kernel_values, mode="full")[
+                        :n
+                    ]
+
+            # Replace NaN/Inf with large but finite values to avoid downstream NaN issues
+            if not np.all(np.isfinite(result)):
+                result = np.nan_to_num(result, nan=1e6, posinf=1e6, neginf=-1e6)
 
             forcing.loc[:, col_name] = result
 
diff --git a/tests/test_modpods.py b/tests/test_modpods.py
index 061c2f8..a0cf482 100644
--- a/tests/test_modpods.py
+++ b/tests/test_modpods.py
@@ -625,8 +625,8 @@ def test_underdamped_kernel_defaults() -> None:
     k = modpods.UnderdampedOscillatorKernel()
     assert k.num_params == 2
     assert k.param_names == ["zeta", "omega_n"]
-    assert k.default_bounds[0, 0] < 0
-    assert k.default_bounds[0, 1] < 1
+    assert k.default_bounds[0, 0] > 0  # Now strictly positive for stable delay dynamics
+    assert k.default_bounds[0, 1] > 1  # Allows overdamped (zeta > 1)
 
 
 def test_kernel_fn_shape() -> None:
@@ -665,18 +665,19 @@ def test_exponential_growth_kernel_defaults() -> None:
     k = modpods.ExponentialGrowthKernel()
     assert k.num_params == 1
     assert k.param_names == ["rate"]
-    assert k.default_bounds[0, 0] > 0
-    assert k.default_bounds[0, 1] > 0
+    assert k.default_bounds[0, 0] < 0  # Negative for stable delay dynamics
+    assert k.default_bounds[0, 1] < 0  # Negative for stable delay dynamics
 
 
 def test_exponential_growth_kernel_increasing() -> None:
-    """Exponential growth kernel should produce monotonically increasing values."""
+    """Exponential growth kernel should produce monotonically decreasing values for stable delay (rate < 0)."""
     t = np.arange(0, 100, 1.0)
     k = modpods.ExponentialGrowthKernel()
-    h = k.kernel_fn(t, 0.5)
+    # Use negative rate for stable delay dynamics
+    h = k.kernel_fn(t, -0.5)
     assert np.all(
-        np.diff(h) > 0
-    ), "exponential growth kernel should be strictly increasing"
+        np.diff(h) < 0
+    ), "exponential growth kernel with negative rate should be strictly decreasing for stable delay"
     assert np.isclose(np.sum(h), 1.0), "kernel should sum to 1"
 
 
@@ -782,7 +783,7 @@ def test_delay_io_train_nse_above_zero(simple_lti_data: pd.DataFrame) -> None:
             windup_timesteps=0,
             init_transforms=1,
             max_transforms=1,
-            max_iter=10,
+            max_iter=20,
             poly_order=1,
             verbose="warnings",
             seed=42,
@@ -1528,3 +1529,155 @@ def test_dual_annealing_reproducible_with_seed(
         model2[1]["kernel_params"].loc[(1, "shape"), :].values, dtype=float
     ).flatten()
     np.testing.assert_allclose(shape_1, shape_2, rtol=1e-10)
+
+
+def test_decoupled_lti_stabilizes_unstable_system() -> None:
+    """decoupled_lti_train should identify an unstable pole without diverging."""
+    np.random.seed(42)
+    n, dt = 500, 0.05
+    T = np.arange(0, n * dt, dt)
+
+    A_true = np.array([[0, 1, 0], [-16, 4.7, 0], [0, 0, -1]])
+    B_true = np.array([[0], [1], [0]])
+    C_true = np.array([[1, 0, 0]])
+    sys_true = ct.ss(A_true, B_true, C_true, 0)
+
+    u = np.zeros((n, 1))
+    u[100:150, 0] = 1.0
+    response = ct.forced_response(sys_true, T, np.transpose(u))
+    u_arr = np.asarray(response.inputs).flatten()
+    y_arr = np.asarray(response.outputs).flatten()
+
+    df = pd.DataFrame(index=T, data={"u": u_arr, "x": y_arr})
+
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        model = modpods.decoupled_lti_train(
+            df,
+            dependent_columns=["x"],
+            independent_columns=["u"],
+            windup_timesteps=20,
+            max_states=3,
+            max_iter=50,
+            verbose="warnings",
+            optimization_method="differential_evolution",
+            seed=42,
+        )
+
+    assert isinstance(model, dict)
+    assert 1 in model
+    result = model[1]
+    assert "A" in result
+    assert "nse" in result
+    assert result["nse"] > 0.5, "NSE {:.4f} is too low for unstable system".format(
+        result["nse"]
+    )
+    eigvals = np.linalg.eigvals(result["A"])
+    assert np.any(np.real(eigvals) > 0), "Expected at least one unstable eigenvalue"
+
+
+def test_observer_compensator_stabilizes_unstable_system() -> None:
+    """decoupled_lti_train + observer-based compensator should stabilize the spring pushcart."""
+    np.random.seed(42)
+    n, dt = 500, 0.05
+    T = np.arange(0, n * dt, dt)
+
+    A_true = np.array([[0, 1, 0], [-16, 4.7, 0], [0, 0, -1]])
+    B_true = np.array([[0], [1], [0]])
+    C_true = np.array([[1, 0, 0]])
+    sys_true = ct.ss(A_true, B_true, C_true, 0)
+
+    u = np.zeros((n, 1))
+    u[100:150, 0] = 1.0
+    response = ct.forced_response(sys_true, T, np.transpose(u))
+    u_arr = np.asarray(response.inputs).flatten()
+    y_arr = np.asarray(response.outputs).flatten()
+
+    df = pd.DataFrame(index=T, data={"u": u_arr, "x": y_arr})
+
+    with warnings.catch_warnings():
+        warnings.simplefilter("ignore")
+        model = modpods.decoupled_lti_train(
+            df,
+            dependent_columns=["x"],
+            independent_columns=["u"],
+            windup_timesteps=20,
+            max_states=3,
+            max_iter=50,
+            verbose="warnings",
+            optimization_method="differential_evolution",
+            seed=42,
+        )
+
+    assert isinstance(model, dict)
+    assert 1 in model
+    result = model[1]
+    A_id = result["A"]
+    B_id = result["B"]
+    C_id = result["C"]
+    D_id = result["D"]
+
+    A_id = np.asarray(A_id, dtype=float)
+    B_id = np.asarray(B_id, dtype=float)
+    C_id = np.asarray(C_id, dtype=float)
+    D_id = np.asarray(D_id, dtype=float)
+
+    # Ensure we have 2D arrays
+    if A_id.ndim == 1:
+        n_states = int(np.sqrt(len(A_id)))
+        A_id = A_id.reshape(n_states, n_states)
+    if B_id.ndim == 1:
+        B_id = B_id.reshape(-1, 1)
+    if C_id.ndim == 1:
+        C_id = C_id.reshape(1, -1)
+    if D_id.ndim == 1:
+        D_id = D_id.reshape(1, -1)
+
+    n_states = A_id.shape[0]
+    assert A_id.shape == (n_states, n_states)
+    assert B_id.shape == (n_states, 1)
+    assert C_id.shape == (1, n_states)
+    assert D_id.shape == (1, 1)
+
+    Q = np.eye(n_states)
+    R = np.eye(1)
+
+    try:
+        K, S, E = ct.lqr(ct.ss(A_id, B_id, C_id, D_id), Q, R)
+        K = np.asarray(K, dtype=float).reshape(1, -1)
+
+        desired_poles = np.array([-4.0, -3.0 + 2.0j, -3.0 - 2.0j])
+        if len(desired_poles) != n_states:
+            desired_poles = np.concatenate(
+                [desired_poles, np.full(n_states - len(desired_poles), -5.0)]
+            )[:n_states]
+
+        try:
+            L = ct.place(C_id.T, A_id.T, desired_poles).T
+            L = np.asarray(L, dtype=float)
+        except Exception:
+            try:
+                _, L, _ = ct.lqe(
+                    ct.ss(A_id, B_id, C_id, 0.01 * np.eye(n_states)),
+                    0.1 * np.eye(n_states),
+                    1e-4 * np.eye(1),
+                )
+                L = np.asarray(L, dtype=float).reshape(-1, 1)
+            except Exception:
+                L = np.zeros((n_states, 1))
+
+        A_cl = np.block(
+            [
+                [A_id - B_id @ K, B_id @ K],
+                [L @ C_id, A_id - B_id @ K - L @ C_id],
+            ]
+        )
+
+        cl_eigvals = np.linalg.eigvals(A_cl)
+        max_real_part = float(np.max(np.real(cl_eigvals)))
+        assert max_real_part < -0.01, (
+            f"Observer-based compensator failed to stabilize the system. "
+            f"Max closed-loop eigenvalue real part: {max_real_part:.4f}"
+        )
+    except Exception as e:
+        raise AssertionError(f"Observer-based compensator design failed: {e}")