Autocatalytic Hysteresis Memory Cell

Use this mechanism as a recurrent cell or continuous-time state-space block. For each hidden coordinate x_i greater than or equal to zero, integrate the input-dependent dynamics dot{x}_i = alpha_i x_i squared plus gamma_0,i plus U_i u_t minus x_i cubed minus beta_i x_i plus W_i h_t plus sqrt(2D_i) times Gaussian noise. Use softplus coordinates or a projected Euler step to enforce nonnegativity, and decode the resulting state together with ordinary features. Pseudocode is: compute gamma_t = gamma_0 + U u_t; compute f_t = alpha elementwise-multiplied by x_t squared plus gamma_t minus x_t cubed minus beta elementwise-multiplied by x_t plus W h_t; sample x_{t+1} = max(0, x_t + Delta-t f_t + sqrt(2D Delta-t) epsilon_t); pass x_{t+1} to the decoder; backpropagate through the unrolled integrator. Initialize alpha, beta, and gamma_0 by numerically finding three roots of alpha x squared plus gamma_0 minus x cubed minus beta x equal to zero and requiring negative drift derivative at the two outer roots. The root calculation comes from the paper; the effective relaxation time tau approximately equal to minus one divided by f_x(x*) is estimated empirically along trajectories. First test synthetic sequence classification with slowly switching regimes and permuted-MNIST, comparing against a GRU and stable linear SSM at equal parameter count. Drive gamma sinusoidally across logarithmically spaced frequencies, record x and gamma, and estimate A with a trapezoidal line integral. The quantitative prediction is a unimodal loop-area curve whose maximum occurs near omega tau equal to one; the measured peak should be within a factor of two of the relaxation-time estimate, while area decreases at both frequency extremes. Increasing D should initially increase switching near the bistable boundary and then destroy coherent hysteresis at large noise, producing a noise turnover.