Generalized Orders of Magnitude for Scalable, Parallel, High-Dynamic-Range Computation
Authors: Franz A. Heinsen
Leo Kozachkov
Venue: Accepted by TMLR ()

Abstract: Many domains, from deep learning to finance, require compounding real numbers over long sequences, often leading to catastrophic numerical underflow or overflow. We introduce generalized orders of magnitude (GOOMs), a principled extension of traditional orders of magnitude that incorporates floating-point numbers as a special case, and which in practice enables stable computation over significantly larger dynamic ranges of real numbers than previously possible. We implement GOOMs, along with an efficient custom parallel prefix scan, to support native execution on parallel hardware such as GPUs. We demonstrate that our implementation of GOOMs outperforms traditional approaches with three representative experiments, all of which were previously considered impractical or impossible, and now become possible and practical: (1) compounding real matrix products {\em far} beyond standard floating-point limits; (2) estimating spectra of Lyapunov exponents in parallel, {\em orders of magnitude faster} than with previous methods, applying a novel selective-resetting method to prevent state colinearity; and (3) capturing long-range dependencies in deep recurrent neural networks with {\em non-diagonal recurrent states, computed in parallel via a prefix scan, without requiring any form of stabilization}. Our results show that our implementation of GOOMs, combined with efficient parallel scanning, offers a scalable and numerically robust alternative to conventional floating-point numbers for high-dynamic-range applications.

Published in Transactions on Machine Learning Research (09/2025)
Generalized Orders of Magnitude for Scalable, Parallel, High-Dynamic-Range Computation
Franz A. Heinsen
frana@ylassvoom.com
GlossRoom Software LLC
Leo Kozachkov
leokoz8 Gagnail com
Thomas J. Watson Research Center, IBM Research
Reviewed on OpenReview: Attps: // openreview. net/ por/formm? id= S0uzb0S00r
Abstract
Many domains, from deep learning to finance, require compounding real numbers over long sequences, often leading to catastrophic numerical underflow or overflow. We introduce generalized orders of magnitude (GOOMs), a principled extension of traditional orders of magnitude that incorporates floating-poin point numbers as a special case, and which in practice enables stable computation over significantly larger dynamic ranges of real numbers than previously possible. We implement GOOMs, along with an efficient custom parallel prefix scan, to support native execution on parallel hardware such GPUs. We demonst that implementation of GOOMs outperforms traditional approaches with three representative experiments, all of which were previously considered impractical impossible, and now become possible and practical: (1) compounding real matrix products far beyond standard floating-point limits; (2) estimating spectra of Lyapumov exponents in parallel, orders of magnitude faster than with previous methods, applying a novel selective-reset method to prevent state colincarity; and (3) capturing long-range dependencies in deep recurrent neural networks with non-diagonal recurrent states, computed in parallel via prefiz scan, without requiring any form of stabilization. Our rosults show that our implementation of GOOMs, combined with efficient parallel scanning, offers scalable and numerically robust alternative to conventional floating-poit point numbers for rygh-dynamic-range applications
1 Introduction
Scientists and engineers often work with real numbers spanning large dynamic ranges, which can exceed the limits of common floating-Por point formats. A typical example is a chain of real-valued matrix products that fails with catastrophic numerical error because it compounds element values beyond reprosentable bounds. Such chains are ubiquitous in science and engineering. For instance, in deep learning, chains of gradients are multiplied together for backpropagation (LeCun et al. 1988). and whether or not these chains explode or vanish determines training success (Hochreiter & Schmidhuber 2022). In control theory, adjoint equations are used to determine optimal control inputs (Bryson 2018). dynamical systerss theory, chains of Jacobian matrix values are iteratively applied to estimale Lyapunov exponents (Strogatz 2018 Pikovsky & Politi 2016). In numerous fields, including economics and finance, Markov models describe behavior over time with chains of stochastic matrices (Ross) 1995..
1997 Pascanu et al., 2013 Mikhaeil et al.
Here, we propose generalixing the concept of "order of magnitude" to include the subset of the complex plane that exponentiates clementwise to the real number line, enabling us to represent any real number- positive, zero, or negative as a complex logarithm that exponentiates to it. We call such complex logarithms "generalized orders of magnitude," or GOOMs for short. As with ordinary orders of magnitude, GOOMs
"The author is Brown University, Providence, RI. 1 Source code replicating our experiments available at githbe..orrkkaasrrre///eeeeeiddd www.ordddrreeo arders_of_magniitude magnitude

---

Published in Transactions on Machine Learning Research (09/2025)
Length of Longest Chain Before Numerical Catastrophe vs. Matrix Size
10a
Longest Chain Length
10a
10
10
10a
101
16
32
64
128
256
512
1024
Square Matrix Size
LMME over GOOMs, Complex64
Conventional MatMul, Flost64
Conventional MatMul, Float32
Figure 1: Longest chain multiplying random normal square matrices without catastrophic numerical error Nvidia GPU, up to maximun of 1M steps, with real numbers represented as Float32 and Float64, and with GOOMs represented as Complex64, applying function we call log- Each point represents the mean of 30 runs, with vertical error bars indicating the standard error across those runs. For every run, at each step, sample matrix elements independently from (0, 1).
or "LMME" (subsection 3.2)
are more stable than the real numbers to which they exponentiate. Formally, we define GOOMs as set of mathematical objects, and show that floating point numbers are a special case of GOOMs. Our work builds upon prior work on logarithmic number systems, with early rools in digital signal processing (Kingsbury & Rayner 1971 Alsuhli et al. 2023 Kouretas & Paliouras 2018 Sanyal et al. 2020y, though the idea of in scientific notation, is older. Defining and naminng GOOMs
representing numbers with exponentiation, enables us to reason about all possible special cases, including floating-poin numbers, in the abstract. From a practical standpoint, GOOMs are complementary to existing numerical formats, not a replacement for them: GOOMs provide a mechanism that can leverage existing numerical formats to enable software applications to operate over a far greater dynamic range of real numbers than previously possible.
We implement GOOMs, well common operations them, including real-valued matrix multiplication, for widely used software framework, PyTorch (Paszke al. 2019), and publish implementation open-source soflware library. Our implemenlation represents real and imaginary components with either Float32 or Float64 numbers, which, as we mention above, are themselves special cases of GOOMs, in effect forming an "edifice of GOOMs." By design, the dynamic range of our implementation far excoeds the dynamic range of Float32 and Floa164, as well as the dynamic range of extended floating-point formats, such as Posits (Gustafson Yonemoto, (017)), proposed as a replacement. for conventional ones (in any practical configuration that could replace floating-poin formats). We also quantify our implementation's relative error, rumning time, and memory use compared to Float32 and Float64. More implementations of GOOMs are evidently possible (including, say, implementations that represent real components with a Posit and imaginary ones with single bit), but we do not explore them here.
We verify that our implementation works as expected with three representative experiments, all which were previously considered impractical or impossible, and now become possible and practical: First, we compute chains of up to 1M random normal square matrices size ranging from 8x8 to 1024x1024, with element values that compound to magnitudes far beyond the limits representable with Float32 or Float64, but which remain representable datatype commonly known as "Complex64." When we compute the chains conventionally over Float32 or Float64, they fail early with catastrophic numerical error, as expected. When we compute the chains over Complex64 GOOMs, they complete successfully (Figure
GOOMs with Float32 real and imaginary components
composite
Second, we formulate, implement, and test parallel algorithm for estimating the spectrum of Lyapunov exponents of dynamical systems, via a prefix scan (Blelloch) 1990), leveraging the greater dynamic range of GOOMs to prevent catastrophic numerical error. We test parallel algorithm on all dynamical systems
2

---

Published in Transactions on Machine Learning Research (09/2025)
from dataset spanning multiple scientific disciplines, including astrophysics, climatology, and biochemistry pin 2023a bl. On parallel hardware, on all systems, our method computes accurate estimates orders of magnitude faster than previous methods, which are sequential (Figure 3b. A key component of our algorithm is a method we devise for conditionally resctting interim states to arbitrary values, as we compute all states in parallel via a prefix scan. We use our selective-reset method to detect whenever interim deviation states are close to collapsing into colinear vectors in the direction of the eigenvector associated with the largest Lyарunov exponent, and to reset such near-colinear deviation states by replacing them with orthonormal vectors in the same subspace, as we compute all deviation states in parallel, over GOOMs.
Finally, we formulate, implement, and successfully train and test on several tasks deep recurrent neural network (RNN) whose layers capture sequential dependencies with non-diagonal state-space model (SSM), computed over GOOMs, executable in parallel via a prefix scan, allowing recurrent state magnitudes to fluctuate freely over time steps, without any form of stabilization. Parallel prefix scans have become more
common in deep learning with the introduction of linear SSMs (Gu ct al. 2021 2022, Smith et al.) 2023 Gu & Dao 2023: Yang et. al. 2024 Feng el. al. 2024 Grazzi el al. 2024). More recent work has extended the use of parallel prefix scans to iteratively compute the trajectories of nonlinear state-space models (Lim et al. 2024 Gonzalez et al. 2024), broadening their potential applicability across more scientific and engineering domains. Our RNN extends the application of parallel prefix scans to non-diagonal SSMs with recurrent states whose elements can fluctuate freely over a greater dynamic range of real values than previously possible, rendering all forms of stabilization unnecessary.
2 Generalized Orders of Magnitude
We define the generalized orders of magnitudes, or GOOMs, as the subset of the complex plane, C' E which exponentiates elementwise to the real number line:
C,
== = (z' C exp(x') RR.
(1)
treating all elements of C' whose exponentiated values are equal to each other as GOOM representations of the same real number. For example, cxp(3 + 2mr) = exp(3 + 4Ti), so we treat 3 + 2mi and 3 + 4ri as GOOM representations of the same real number, exp(3) = 20.0855.
Equivalently, we can say that GOOM's imaginary component must be either a non-zero integer multiple of mi integer multiple of 2mi, such that the imaginary component's exponentiation is in {-1,1), 1), per Euler's identity. When the imaginary component is an integer multiple of 2Ti, the real number to which the GOOM exponentiates is positive; otherwise, the imaginary component is a non-zero integer multiple of Ti, and the real number represented by the GOOM is nogative. As a convention, we treat zero in the real number line as non-negative-icc. we represent it as a positive GOOM. Ordinary orders of magnitude, including tog-probabllities, are in the subset of GOOMs with zero imaginary component.
The concept of complex order of magnitude, defined as x' log for x E I, may initially appear unconventional; however, it is well-defined since ex x. We stop short of formally inducing isomorphism from IR to C², and vice versa, because we want to keep open the option of applying arbitrary complex-values transformations that leverage the structure of the complex plane 2 Our only requirement for a number to be a GOOM in C' is that it must exponentiale to a number in the real line.
As with ordinary orders of magnitude, GOOMs are more stable than the real numbers to which they exponen- tiate. We now provide two examples of how GOOMs can be used to perform two standard operations- scalar multiplication and the dot product of two real-valued vectors in way that is inherently more numerically stable than traditional methods.
"Consider, for example, deep learning model that processes data in C and includes final layer that applics a transformation from C to C' thereby allowing the data be scaled and projected R. Not coincidentally, our implementation of GOOMs makes instantiating and training such model ttraightforward, because it ensures that backpropagation works seamlessly over C, over C, and across mappings between C and RR, by extending Poror'''s 's infrastructure for complex data types.
3

---

Published in Transactions on Machine Learning Research (09/2025)
Example 1 (Scalar Multiplication in R Becomes Addition in C²) Suppose given d real num- bers 22:22: .. Ld, and we wish to compute. their product:
d
y = x1 22
Dd
zyy.
j=1
We can express this product as a sum in C' since, denoting x'J = log 2y for j = 1,2,...,d
exp (1)) II exp (2y)
y
Example 2 (Dot Product in IR Becomes log-sum-exp in C¹) Suppose we are given two d-dimensional vectors and with elements in R. We wish compute their dot product:
c= a
ag bj.
We can express this dot product as log-sum-eap (LSE) operation in C' Denoting e' as the element-wise sum of a' = log(()) and 8' = log(b),
= + b/j,
we find that d' log(c)== LSE LSE (z'). This is because
exp(x') Σ exp(a,) exp(t)) Σ exp(xs)
logoge log Σ exp(2)) LSE(2').
1=1
Note that even if the individual elements of and b very large, the GOOM representation in C' remains numerically stable. For example, suppose aj = exp(1000) Naively computing the dot product in R fails due to floattin-poois point overflow. However, computing the LSE in C works easily, since 2/ 2000, which is well within the capabilities of the numerically stabilized implementation of this operation, e.g., in PyTorch.
2.1 Floating-Poii Point Numbers are Special Case of GOOMs
Given a GOOM x' = a + bi, representing real number x, we know from equation 1 that x = e² eatbi ea x bi E R. Calling e the "base," the "exponent," and bi which is in {-1,1} 1}, the "sign," we have:
exponent
x E R.
(2)
base
sign
If equation 2 seems familiar, it is because floating point representations of real numbers, which have a base, an exponent, and a sign, are in C² Their base is 2 instead of e, their exponent, is encoded in two sequences of bits (one representing a significand or mantissa, the other reprosenting an integer exponent) instead of as the real component a complex number, and their sign encoded as a bit instead as an exponentiated imaginary component, but otherwise they represent real numbers the same domain. Given floating point number with signed significand s, represented by fixed number of bits, and integer exponent n, also represented by fixed number of bits, we have:
œ x 2n
abs(s) x 2d x sign(s)
(3)
2108, (abs(e)) 2n x sign(s)
(abs(()) x sign(s), x E IR.
4

---

Published in Transactions on Machine Learning Research (09/2025)
In other words, floating-point numbers are a special case of GOOMs Unlike floating-poin numbers, which are defined in terms of a fixed number of bits, GOOMs are defined more generally in terms of real and imaginary components, independently how those components might be represented in a computer. Also unlike foating-point numbers, GOOMs explicitly leave open the possibility of applying arbitrary transfor- mations that leverage the structure of the complex plane. In practice, GOOMs are complementary to existing floating-point formats, providing a mechanism that enables software applications to operate over a greater dynamic range of real numbers than previously possible, without requiring the introduction of and support for new moating-point formats, e.g., in hardware devices like Nvidia GPUs.
3 Implementation
We implement GOOMs for PyTorch, widely used software framework for parallel computation, as Com- plex64 and Complex128 composite data types, which represent their real and imaginary components Float32 and Float64 numbers, respectively. That is, we leverage existing foating-poin formats, which are special cases of GOOMs (2.1), to represent the real and imaginary components of other GOOMs, for rep- resenting real numbers over a greater dynamic range than is possible with existing floating-point formats. In practice, implementing this "edifice of GOOMs" is more straightforward than our description suggests, thanks to the availability of and support for complex data types in PyTorch.
The dynamic range of our implementation is, expected, far greater than that of Float32 and Float64 (Table [), the two numerical formats with highest dynamic range supported by Nvidia GPUs. It is also far greater than that of extended floating-point formats such as Posits in any practical configuration that could replace Flooa322 or Float64 | We also note that our implementation can be readily extended to represent real and imaginary components with Posits or other numerical formats that make different tradcoffs compared to Float32 and Float64, provided such formats are supported by hardware. Other implementations are evidently possible (including, for example, those that represent real components with a Posit and imaginary ones with a single additional bit) but we do not explore them here. It bears repeating that in practice GOOMs are complementary to existing numerical formats.
Representation
Bits
Smallest Normal Magnitude
Largest Normal Magnitude
Float32
32
10-38 38 = exp (- 13399333
10
1038 = exp (101.9395)
Float64
64
10-308 308 ~ exp (800
10008 ~ exp (102.8506)
Complex64 GOOM
64
exp ( - 1038)
exp (1038)
Complex128 GOOM
128
exp (- 303
exp (10308)
Table 1: Dynamic range for Complex61 and Complex128 GOOMs versus Float32 and Float64.
By design, Complex64 and Complex128 GOOMs benefit from the greater precision at smaller magnitudes of Float32 and Float64, respectively, for representing larger magnitudes. Float32 and Float64 represent magnitudes between 0 and 1 with negative exponents, consuming approximately half of all possible exponents in the floating-point format, and magnitudes above 1 with positive exponents, consuming the remaining exponents, such that precision decays as magnitude increases. Complex64 and Complex128 GOOMs thus have greater precision at smaller magnitudes (in fact, the magnitudes become too small to be mapped via exponentiation to Float32 and Float64, respectively), and precision decays below that of the underlying point format as magnitude increases toward the format's maximum representable magnitude, and beyond (Figure 2).. We compare the relative error, running time, and memory use of Complex64 and Complex128 GOOMs versus Float32 and Float64 in Appendix D
a special case too. 4 For example, the 64-bit confiurration of Posits proposed as a drop-in replacement for Float64 by Gustafson & Yonemoto 2017) (Table 3) covers magnitudes from roughly 10-299 exp(-102.8576 10299 exp(102.8336, for less than our imple- mentation of Complex64 GOOMs, which cover from roughly exp(-1038) to exp(1038) excluding subnormal real components.
Extended floating- point numbers, including Posits (uustafson & Yonemoto 2017
5

---

Published in Transactions on Machine Learning Research (09/2025)
0
1
e
C
55
15
55
+>
C
0 1
log c
C
1
Figure 2: At the top, show the range of magnitudes with positive sign representable by Float32 Float64, up to a maximum c, and their approximate share n, the number of possible bit sequences. For Float32, = 2" for Float64, 2 At the boltom, show the represented by the same ploating-pinll format. The shares of n are approximate to account for bitwise differences between Float32 and Float64. For magnitudes with negative sign, the diagram identical.
magnitudes, mapped
complex GOOM's real component,
Along with functions for mapping existing floating-poins numbers to complex-typed GOOMs, and vice versa, we also implement variety of real-valued functions commonly applied in scientific and engineering applications, including matrix multiplication, over complex-typed GOOMs. As expected, such functions are more stable over C. For example, scalar product over R becomes scalar addition over C (Example I sum over R bocomes a log-sum of exponentials over C² i.e., it remains a sum, but with elementwise nonlinear transformation ions before and after the sum (Example 2 A sum of products over R, as a matrix product or tensor contraction, bocomes a log-sum-cxp of scalar additions over C', i.e., a sum of sums, with elementwise nonlinear transformations before and after the first sum. A composition of linear transformations and elementwise activation functions, in deep neural network, becomes sum of sums of sums, with elementwise nonlinear transformations interspersed between sums.
All implemented functions are parallelized, broadcastable over an arbitrary numbers of preceding, indices, and compatible with backpropagation of gradients over C, over C', and over mappings between and R, taking special care to handle the singularity at zero gracefully, for use in a broad range of applications, including deep learning.
3.1 Mapping Between Floating-Poid PPint Numbers and GOOMs
Floating-Point Numbers to Complex-Type Typed GOOMs Given a floating-point number x, we map it complex-typed GOOM s' as follows:
x' -- log(x:)
real
0,
x 2 0 < 0. ni.
(4)
log(r) = log
abs (z:)
+
x
real where abs(-) and log(-) denote custom implementations of real absolute value and elementwise logarithm, rospectively, necessary to accommodate broad range of applications.
real
rel Our implementation of abs(-) redefines its elementwise finite derivative, with respect its input, always be non-zero, treating input values equal to zero non-negative by convention,

---

Published in Transactions on Machine Learning Research (09/2025)
d (r)
x20 0 -1, 0,
1,
redefined finite derivative
(5)
dx
ensuring that gradients are zero in backpropagation only when the backpropagated error is zero. Otherwise, our implementation of real absolute value behaves like PyTorch's default implementation.
Our implementation of log(-) is configurable, such that for real input elements equal or numerically close to zero, we can specify whether the function will (a) return a sentinel value representing negative infinity, maximizing precision up the limits of the chosen data type; or (b) return finite floor value that numerically exponentiates to zero, for applications that require values always to be finite
roul We redefine the elementwise finite derivative of log(:), with respect to its input, to add data-type-specife small number to its denominator,
d llg(2))
If redefined finite derivative
(6)
dx
x t
ensuring that gradients always finite for backpropagation, including at the singularity for the logarithm of zero. Otherwise, our implementation of real logarithm behaves like Py default implementation (e.g., it is undefined for negative real input values).
Complex-Typed GOOMs to Floating-Pore Point Numbers Given complex-typed GOOM x', we map it to a floating-point number æ as follows:
x + exp(x')
exp(x') = R exp (x'))
where exp (-) denotes a custom implementation of complex elementwise exponentiation, and R(-) denotes real component. We discard the imaginary component, which may not always be zcro duc to accumulation of numerical errors. Recall that GOOMs exponentiate to R, by definition.
We redefine the clementwise finite derivative of comples (.), with respect to its input, to always be non-zero for backpropagation, by adding to or subtracting a data-lype-specific small number that shifts the finite derivative's real component away from zero:
complex CXP (x')
E,
R
complex
EXP (x')
2 0
comptex
exp
(x') +
// redefined finite derivative
dr
C,
R
complex CXp
(x')
<0
ensuring that gradients zero in backpropagation only when the backpropagated error is zero. Otherwise, our implementation of complex exponentiation behaves like PyTor'''s 's built-in implementation, relying on Py default behavior for backpropagating gradients over complex data types.
3.2 Real-Valued Matrix Multiplication over Complex-Typed GOOMs
Given two matrices, A' Cm and B' Cddxm the equivalent of their matrix product over IR is expressible log-sum-exp of elementwise additions:
IIn initial implementation, specify (og(SNN)) the fiuite floor value, where SNN is the smallest normal number representable by the chosen floating-point format. For example, for Complex64 GOOMs, which have Float32 components, the smallest normal number representahle by Float32 approximately 1.18 10-38, so the finite floor, log(1.18 10 10-38) -36) is approximately -174.7, which numerically exponentiates to at Float32 precision.

---

Published in Transactions on Machine Learning Research (09/2025)
LMME (A',B) = log exp(A') exp(B')
MatMul over
log Σ exp(Alg) & exp(
(9)
R xmx
LSE J
Aly B's
Crexdx
where LMME is shorthand for "ogg-matri--mlltthieetioo-appppp LSE denotes log-sum-cxp, and @ and @ denote elementwise product and addition, respectively, of each row of the first matrix with each column of the second one (i.e., X and + denote outer product and addition, respectively, over index j). The final expression in equation g is casily parallelizable, bocause each row-column pair can be clementwise summed- tten--o--sum-expoeenttafrt independently of all other pairs.
An optimal parallel implementation of LMME requires an efficient parallel implementation of log-sum-exp of pairwise elementwise sums. A straightforward approach to implement it would be to compute all pairwise elementwise sums in parallel, broadcasting over all other dimensions, then apply log-sum-exp, but doing would be impractical, because it would require O(ndm) space. Another obvious approach would be to apply log-sum-exp to the elementwise addition of each pair of vectors independently of the other pairs (e.g., with a poctor-mapping, or "vmap," " operator), but doing so would run into menory-bandwidth constraints on hardware accelerators like Nvidia GPUs, which are better suited for parallelizing computational kernels that execute and aggregate results over tiled sub-tensors. Unfortunately, Py and its ecosystem, including intermediate compilers like Triton (Tillet et al. 2019) currently provide no support for developing highly optimized complex-typed typed kernels.
As compromise, we implement LMME(A', B') so it delegates the bulk of parallel computation to PyTorch's existing, highly optimized, low-level implementation of the dot-product over R, as follows:
Compromise
LMME (A', B') log exp (A/ - a.) exp (B', bE) + as + bb.
(10)
Scaled MatMul
where E IR" and by Rh are vectors with log-scaling constants, necessary because the interim exponen- tiation to R may return values outside the bounds representable as floating-point numbers:
elem max
max (R (A(y)), 0)
=
j
(11)
elem 6% = max
max (RR (B)x)) 0).
with m-x denoting elementwise maximum, and R(-) denoting elementwise real components. The expression in equation I0 follows from the fact that we scale each row of the left matrix and each column of the right matrix, before multiplying the two matrices, and subscquently undo the scaling for every element of the resulting product:
elem
eepp (B'x)) )) exp(ba)
exp (A(y) exp(a:)
Σ exp (A',) exp (B'r) =
exp
Σ
c exp(a;) O exp(bx),
(12)
where division applied clementwise, implicitly broadcasting over missing indices, and @ denotes elementwise (Hadamard) product, also implicitly broadcasting over missing indices. We compute a, and by without impacting gradients in backpropagation (f.e,, detached from the computation graph).

---

Published in Transactions on Machine Learning Research (09/2025)
We recognize that our initial implementation of LMME is sub-optimal compromise, both in terms of precision (we execute scaled dot-products over a floating-point format, instead of elementwise sums over complex-typed GOOMs) and performance (we must compute not only a scaled matrix product, but also per- row and per-column maximums on the left and right matrices, respectively, two clementwise subtractions, and two elementwise sums). In practice, we find that our initial implementation of LMME works well in diverse experiments (section D. incurring exccution times that are approximately twice as long as the underlying real-valued matrix product on highly parallel hardware a reasonable initial tradeoff, in our view, for applications that must be able to handle greater dynamic range of real magnitudes. See Appendix D for comparison of our implementation of LMME to PyTorch's s matrix product over Float32 and Float64.
3.3 Other Real-Valuee Functions over Complex- Typed GOOMs
In principle, can naively lormulate the equivalent over C' of any real-valued function ʃ as function composition f that elementwise exponentiates the complex input, applies f to the exponentiated input, and then takes the elementwise logarithm of f's output:
I' = log o f exp.
(13)
In practice, we can never implement the naive formulation, because the interim exponentiation to R may return values that fall outside the bounds representable by existing floating-point formats. Instead, we must either (a) implement f' such that it avoids interim exponentiation to altogether, always remaining in C', or (b) scale elements in the log-domain, before exponentiation to E, and undo the scaling after taking the elementwise logarithm, using techniques analogous to the "log-sum-exp trick," " but specific to each function f. We use both approaches, (a) and (b), to implement over C' a variety of real-valued functions commonly used in science and engineering.
Please see our published source code for the full list of implemented functions.
4 Representative Experiments
We verify that our implementation works as expected with three representative experiments: (a) long chains of real-valued matrix products that compound element magnitudes far beyond the limits of Float32 and Float64 numbers; (b) parallel algorithm for estimating the full spectra and largest of Lyapunov exponents of dynamical systems, via a prefix scan, orders of magnitude faster than previous methods, leveraging complex-typed GOOMs to avoid catastrophic numerical error; and (c) a deep RNN whose layers capture soquential dependencies with non-diagonal lincar SSM over complex-typed GOOMs, executable in parallel via a prefix scan, allowing recurrent stale magnitudes to fluctuate freely over a greater dynamic range of real numbers than previously possible, making all forms of stabililization unnecessary.
4.1 Stability and Dynamic Range of Complex GOOM Matrix Products vs. Conventional Methods
We compare the stability and dynamic range of our implementation of LMME against those of conventional matrix multiplication over E with floating-poing point numbers, on long chains of matrix products that compound clement magnitudes toward infinity. For in {8, 16.32.... 1024), for t in {1, 2...,.....100 10%), we multiply a chain of random square matrices A, E Rdxd each element of which is independently sampled from a normal distribution (0, 1). At each step t, the state of chain is
S, = A.SL11 A, ~ NO 1001d
(14)
with So N/00, 1)dxd as the initial value. We We update the chain's state until all steps are completed, or until the computation fails with catastrophic numerical whichevre occurs first. For each value of d, attempt to complete the chain 30 times recent. Nvidia GPU.
9

---

Published in Transactions on Machine Learning Research (09/2025)
When we attempt to compute all chains with conventional matrix multiplication over R with Float32 and Float64 numbers, we find that all chains fail early with catastrophic numerical error.
When we attempt to compute the chains over C' with our implementation of LMME, representing GOOMs as Complex64 numbers, we find that all chains successfully complete all steps (Figure !
At each step i, the log-state of chain over C' is,
S LMME (A),),-1), Ai log 10,,
(15)
with S log,V(0, 1)axd as the initial value. As before, we repeatedly update the chain's state until all steps are completed, until the computation fails with catastrophic numerical error- whicheeer occurs first. For each value of d, attempt to complete the chain 30 times recent Nvidia GPU.
We repeat all tests on a recent multi-core CPU, and obtain essentially the same results. We do not show plots for the CPU lests because they are indistinguishable from those executed on the GPU. Further analysis reveals that all chain states are representable their entirety as complex-typed GOOMs, but not as the floating-point numbers to which they would exponentiate, were it possible.
4.2 Parallel Estimation of Lyapunov Exponents over Complex-Typed GOOMs
A fundamental concept the theory of dynamical systems is the Lyapunov exponent (Strogatz) 2018
Pikovsky Politi 2016 Bradley & Kantz, 2015; Gilpin 2024). Loosely speaking, Lyapunov exponents (LEs) quantify how nearby trajectories in the state space of nonlinear dynamical system diverge (or converge) over time. In general, Lyapunov exponents (LEs) must be computed numerically along specific trajectories of the system under considcration. If the system is stable-for instance, if it is contractive (Lohmiller &
Slotine 1998 Srinivasan & Slotine, 2023 it may be possible to conclude a priori that the L.ss are negalive. However, assigning precise numerical values to the LEs cannot in general be done analytically. Instead, LEs must be computed by repeatedly multiplying Jacobian matrices of the dynamics. However, this approach can be numerically unstable if the underlying dynamics themselves are unstable. Methods that address this instability, such iterative QR decomposition (Pikovsky & Politi, 2016), are not parallelizable in time. To understand why, it helps to recall the standard algorithm for computing the largest Lyapunov exponent (cf. (Pikovsky & Politi, 2016, Section 3.1)). The method propagates a vector through a sequence of matrices, normalizing this vector at each step to keep it on the unit sphere. This normalization prevents numerical overflow or underflow. However, bocause the normalization depends on the current state (specifically, its magnilude), the procedure cannot be implemented using a parallel scan (Blelloch 1990).
In particular, consider the nonlinear dynamical system
2, = f((xx-1) E RR,,
(16)
where Je denotes a given sequence of nonlinear functions. A simple example of equation 16 is a recurrent neural network time t. If we wish to understand how small spatial perturbations to solutions of equation 16 influence the future behavior of the system, we may analyze the associated variational equation
tanh(Wz,, + 6 + uz), where W is a weight matrix, b is bias, and 2n is an input at
Of.
8x, Jg d
where
J 0x
(17)
By unrolling the variational equation equation 17 in time, we see that the dynamics of spatial perturbation to the solution of the nonlinear dynamics equation 16 at time t = 0 is entirely determined by the product of Jacobian matrices up to time L:
dr. Jr 520
[[ R=1
Jk
5mm
H, Exr.
where He defined the cumulative product of Jacobian matrices up to time t. Using H,, the LEs are defined the long-time average exponential growth (or decay) rates of these infinitesimal perturbations.
10

---

Published in Transactions on Machine Learning Research (09/2025)
Formally, the Lyapunov exponents X given by
x lim 1000
log
H1)
(18)
where (.) denotes the i-th singular value of the product of Jacobians and log denotes natural logarithm. In other words, they quantify the average exponential rates which volumes spanned by perturbation vectors expand contract under the system's 's dynamics.
10a
Sequential Estimation Time
10a
8 Multiple Paralle
10
10c
Mean, All Systems Dev. All Systems
Number Steps
Figure 3: Time to estimate the spectrum of Lyapunov exponents sequentially, as a multiple of time to estimate it in parallel, increase the number of steps, for all dynamical systems dataset spanning multiple scientific disciplines Gilpin 2023a b). in single Nvidia GPU. The improvement starts tapering off at 10 steps because the GPU's compute capacity saturated by parallel QR decompositions all steps. Appendix A shows plots by system.
4.2.1 Full Spectra of Lyapunov Exponents in Parallel
A1 A2
The standard approach for estimating the full spectra of Lyapunov exponents,
i8:
rst
A
log (abs (Diag ((R,)))
(19)
ATT
where
Q₁, QR(S.)
QR-decomposition
(20)
S. = J.Q.-1
given initial deviation states with unit norms So Rddd That is, al each step / we obtain a triangular matrix R, by QR-decomposition of St, obtained by applying J, to the preceding orthonormal basis Q obtained by previous QR-decomposition. The clementwise logarithms of R, 's eigenvalues, in its diagonal, are unscaled estimates of A at step I; their scaled mean is the estimate of A.
Parallel estimation A requires not only that we handle magnitudes that may not be representable as floating-point numbers, which we already know we can do with complex-typed GOOMs, but also that we execute QR-ecooposition before and after the application of each Jacobian matrix Jr, to obtain its input state's orthonormal basis Q. and its output state's triangular factor RG, respectively. Alas, if we
11

---

Published in Transactions on Machine Learning Research (09/2025)
naively attempt to compute all deviation state matrices in parallel, via a prefix scan, with the intention of subsequently executing all QR decompositions ions in parallel, we are not able to do so, because states tend to collapse i.e., become colinear in the direction of the eigenvector associated with the largest Lyapunov exponent. We must find a way to prevent vectors from collapsing into colincar states.
Our solution an algorithm that can compute all deviation state matrices in parallel, by applying a method we devise for conditionally resetting interim states arbitrary values in any linear recurrence, as we compute it in parallel via prefix scan. Our conditional-resetiin enables to detect whenever any interim deviation states are close to collapsing into colinear vectors, and to resct such near-colincar states by replacing them with orthonormal vectors in the same subspace, at any time step, as we compute all deviation states in parallel, via a prefix scan, over complex-typed GOOMs. We describe our parallel algorithm below.
method, which describe separately in section [
Parallel Algorithm We execute the following groups of parallelized computations sequentially:
(a) Compule input states, S1, S2, ..... Sr-1: I. by cumulatively applying all Jacobian matrices except the last one, J1, J2..... : 1r to the given initial state S., in parallel, via prefix scan, over complex- typed GOOMs (to be able to handle all magnitudes), applying our selective-resetting method (section 5) to reset any interim states in the chain that are close to becoming colinear, replacing them with an orthonormal basis obtained via QR-decomposition which we apply in parallel only to those interim states, if any, selected for resetting. We define "close to becoming colinear" the cosine similarity of any pair of state vectors exceeding a specified threshold. Before applying QR-decomposition to any near-colinear interim states, first log-scale them to log-unit norms, over complex- typed GOOMs, and then exponentiate them to values representable as floating-point numbers.
(b) Compute orthonormal input bases Qo, Q1, Q2 .... Qr-1, 11 by applying QR-decomposition to every first log-scale its input vectors to log-unit
Sy- Before applying each such QR-decomposition, norms, over complex-typed GOOMs, and then exponentiate them to values representable as floating- point numbers. We execute every QR-decomposition independently of the others, in parallel.
(c) Compute output states, Si, Sÿ, ...; Sj, by applying each Jacobian matrix Jr to its preceding input basis Q.-1- We apply each Jacobian matrix to its input basis independently of the others, in parallel.
(d) Compute the estimated A by applying QR-decomposition to every output state Sü, extracting the diagonal elements (eigenvalues) of each triangular matrix R,, taking the elementwise logarithm of their absolute values, and computing their mean over all states. We apply QR-decomposition and extract and transform diagonal elements for each state independently of the others, in parallel.
With T-way parallelism, fixing d, our algorithm's time complexity is (loggT); because the parallel prefix scan in group (a) above has (log I') time complexity, and every subsequent group of computations has O(1) time complexity, except for the final mean computed in group (d), which has Ollog T) time complexity. We test our parallel algorithm all dynamical systems from dataset spanning multiple disciplines, including astrophysics, climatology, and biochemistry (Gilpin 2023a bb, on recent Nvidia GPU, We find that our parallel algorithm computes accurate estimates, orders of magnitude faster than sequential estimation with previous methods (Figure
4.2.2 Largest Lyapunov Exponent in Parallel
Estimating only the largest Lyapunov exponent (LLE) in parallel is straightforward over complex-typed GOOMs, thanks to their ability to handle a greater dynamic range of magnitudes. The standard approach for estimating the LLE relies on measuring changes in norm of an initial deviation vector E Rd with unit norm:
est LLE
T
|s Il
log
(21)
AIT
|
12

---

Published in Transactions on Machine Learning Research (09/2025)
where At is the discrete time interval, and
8t J
(22)
|ssll|
Normalization of each preceding deviation state necessary to keep element magnitudes from compounding above or below the bounds representable by existing floating-poing poini formats. Without normalization of preceding steps, equation 21 becomes:
est
|| ) |
LLE
log
81 J181-1,
80 uo.
(23)
ATT
In Appendix B
show that equation 23
simplifies to:
LLE
LSE(2 (U ub))
(24)
2ATT
where PSCAN(LMME)( denotes a parallel prefix scan applying LMME over time steps, and
For each 2 {uo, J1, J2, ..... JTh, d'
We test our method for parallel LLE estimation on all dynamical systems from dataset spanning multiple disciplines, including astrophysics, climatology, and biochemistry (Gilpin 2023a on recent Nvidia GPU, We find that parallel method computes accurate estimates, and as expected, its execution is orders of magnilude faster than sequential estimation with previous methods.
4.3 Parallelizable Non-diagonal State-Space RNN over Complex- Typed GOOMs
Deep learning models apply variety of normalization methods (RMSNorm, LayerNorm, BatchNorm, Group- Norm, etc.), residual mechanisms (LSTM cells, GRU cells, residual layers, skip connections, etc.), gated lincar activation functions (ReLU, SiLU, GeLU, GLU, SwiGLU, etc.), and various other techniques to keep the magnitudes of feedforward elements and backpropagated gradients from becoming too small or too large, for preventing catastrophic numerical errors. Deep RNN models tend to be most susceptible to such numer- ical issucs, because their layers capture sequential dependencics by repeatedly applying the same recurrent transformation at each lime step, compounding the transformation.
We formulate and implement a deep RNN whose layers capture sequential dependencies with non-diagonal recurrences over complex-typed GOOMs, allowing recurrent state magnitudes to fluctuate freely over greater dynamic range of real numbers than previously possible, enabling compulation in parallel via prefir scan without requiring any form of stabilization. Otherwise, the RNN operates conventionally over floating-point numbers. The RNN consists of a standard embedding layer that maps token values (e.g., integers representing symbols in a vocabulary) to feature vectors representing token states; multiple residual recurrent layers that capture sequential dependencies over complex-typed GOOMs, and conventional model head that task-specific (e.g., linear transformation for classification tasks).
Each residual recurrent layer captures sequential dependencies multiple heads per token, by applying: first, LayerNorm and a linear transformation with bias to obtain each token's 'ike'ss heads; second, a parallel prefix scan of non-diagonal linear recurrence, capturing per-head sequential dependencies with a standard state- space model, (over complex-typed GOOMs, as we will detail shortly); and finally, gated linear units (GLU) on every head, followed by a linear transformation of the flattened heads, to obtain residual, added to the token's 's input state.
Per head, capture sequential dependencies with standard linear state-space model:
13

---

Published in Transactions on Machine Learning Research (09/2025)
Ft Ar1-1 Buu
(25)
Yt = Cr. + Dut.
where E Rd, E IRd, and E RR, with h = 2d for application of GLU, and A, B, C, and D are matrix parameters but we compute the recurrence over compler-ypee GOOMs. That is, we first map the initial state To, all input head states ut, and all matrix parameters to camplex-type GOOMs,
For each {ro, We: A, B,C, C, DJ,
2'
then compute the recurrent relationship in equation 25 over complex- typed GOOMs, for every head, in parallel, via a prefix scan that obtains all states x' without applying any form of stabilization,
x/- LSE x1/
u/)).
(26)
then map of back to floating-point format, for executing the remaining computations in the layer. We cannot exponentiate the elements of z outright, because their magnitudes might not be representable by a floating-poin format, so we first log-scale 5y, then exponentiate its elements, as follows:
e ef max ()(xz))) scaled 24 t- exp (1:: c 2),
(27)
computing the log-scaling real constant without altering gradients in backpropagation (f.e. detached from the computation graph). The elements of every scaled z, fall between cxp(2) and exp(2).
We train and test instances of our RNN several tasks, including (but not limited to) generative language modeling on The Pile [Gao et al. 2020), classification and generation of pixel sequences from MNIST (LeCnn et al. 2010), and a Copy Memorr lask. Perhaps the most remarkable finding about the training dynamics is how unremarkable they are, even though we are computing non- non-diagonal recurrences in parallel without ony form of stabilization. We are even able to compile and autocast Float16 all components of the RNN that execute operations over floating-poin formats, executing in PyTorrch's 's eager mode only the non-diagonal recurrences in parallel, via the prefix scan over complex-typed GOOMs. Figure shows examples of training runs for two tasks. Please see our published source code to replicate training runs on all tasks.
5 Selective-Resettin Method for Parallel Scans of Linear Recurrences
As part of our work to formulate algorithm for parallel estimation of Lyapunov exponents (subsection 4.2 we devise a method for conditionally resetting interim states to arbitrary values in a linear recurrence, as we compute all of its states in parallel, via prefix sean. For ease of exposition, we describe our method with a non-diagonal, timo-variant, linear recurrence over R, understanding that the method generalizes to any linear recurtence diagonal or not, time-variant or not over IR or other fields.
Given:
an initial state Xo in idddd
sequence of transition matrices A,, Ag...... ................................ An, each in pdxd
a "selection" function S Rdxd + {0, 1), and
"reset" function R: Rdxd + Idxd
we wish compute a linear recurrence X₂ A,X,-1, for t in (1,2,..,,n), via a parallel prefix sean, but modifying the recurrence, such that, for any interim compound state A computed during the parallel scan,
14

---

Published in Transactions on Machine Learning Research (09/2025)
Natural Language Generation, Sequences of 1024 Tokens, The Pile
Classification from Last Pixel Value, Sequences of 784 Pixel Values, MNIST
2.5
10
Crozz Entropy Lass
Cross-Entropy Loss
1.5
1.0
0.5
0.0
107
10
10
Number of Tokens
Epoch
Pigure 4: Examples of training dynamics for the RNN implement, capturing sequential dependencies with non- diagonal recurrences, computed in parallel via prefix
without any form of stabilization. Left: Natural language generation on The Pile (Gao al., 2020), with 124M-parameter RNN incorporating a 50257 token-id vocabulary and 24 layers; we stopped training at TOB tokens. Right: Classification, from last pixel value, of sequences of 784 token-id vocabulary layers. See our source code for replicating all training runs, including for both tasks shown here.
pixels from MNIST (LeCun al. 2010), with a 12.8M-parameter RNN incorporating
if ((A)) = 1, we reset A by replacing it with R(A"), making R(A") the new initial state for subsequent states, until we reach either (a) the final state of the recurrence, or (b) subsequent state that has previously been selectively reset whichever occurs first.
We accomplish our objective in two steps. First, we add bias matrix each step: X; = A,X,11 Br, with every bias matrix initialized to all-zeros: B, - [0]axd for & in (1, 2....,,., n), such that, at initialization, the new recurrence matches the original unmodified one.
Second, we execute a parallel prefix scan over the new recurrence, but instead of applying an ordinary linear recurence with bias, we apply the following binary associative transformation to every pair of tuples containing preceding interim states, (Aprev: Bie Bprev), and subsequent ones to be updated, (Acurr Beur)):
1/ first, selective resetting:
If (Aprev) = 1 and Bprev {0]axd
B prev -- R (Apres
Aprev -- {0]axe
(28)
// then, ordinary recurrence:
A curr
A A
B A B Birev B
We obtain modified sequence of states that may may not match the original sequence, because more of its states may have been reset during the scan. The seloctiveresetting transformation equation 28 is associative, because it can reset each state only once, and the zeroed-out transition matrix of any reset state eventually zeroes-out all subsequent transition matrices via cumulative multiplication.
Our selective-reecttnng method incorporates multiple moving parts, and it may take some effort to grasp fully how all of them interact upon a first read. For the convenience of readers seeking an intuitive understanding of the method, Appendix C explains it more informally with step-by-step examples.
15

---

Published in Transactions on Machine Learning Research (09/2025)
6 Discussion
By encoding real numbers as complex logarithms, GOOMs enable stable, scalable and parallel computation across dynamic ranges far beyond what Float32 or Float64 support, without requiring changes to hardware or floating-point standards.
Our experiments support GOOMs' robustness and versatility. long matrix product chains, complex-typed GOOMs avoid overflow and underflow where conventional approaches fail early. In Lyapunov exponent estimation, the combination of complex-typed GOOMs with parallel prefix rosetting method enables accurate, timo-parallel analysis of chaotic systems, orders of magnitude faster than with previous methods. In deep learning, RNNs handle freely fluctuating non-diagonal recurrent states over complex-typed GOOMs without compounding beyond representable limits or degrading gradients, rendering all forms stabilization unnccessary.
and novel selective-
We show that GOOMs can be implemented straightforwardly native PyTorch, supporting aulograd and GPU execution through standard complex types. We note that although the current implementation of operations like mammmmi-mmliippcattooooo (Section 3.2 introduces some overhead, future work custom kernels could close this gap.
To summarize, GOOMs in practice offer robust, software-level solution to numerical instability in high- dynamic-range computations. Their blend of flexibility, precision, scalability, and parallelizability makes them a powerful tool for scientific computing in general, and deep learning particular.
References
Chada Alsuhli, Vasileios Sakellarion, Hani Saleh, Mahmoud Al-Qutayri, Baker Mohammad, and Thanos Stouraitis. Number systems for deep neural network architectures: A survey. arXiv preprint arXiu:2307.00085, 2023.
Guy E. Blelloch.
Prefix sumn and their applications. Technical Report CMU-CS-90-190, School of Computer Science, Carnegie Mellon University, November 1990.
Elizabeth Bradley and Holger Kantz. Nonlinear time-series analysis revisited. Chaos: An Interdisciplinary Journal of Nonlinear Science, 25(9), 2015.
Arthur Farl Bryson. Applied optimal control: optimization, estimation and control. Routledge, 2018.
Leo Feng, Frederick Tung, Mohamed Osama Ahmed, Yoshua Bengio, and Hossein Hajimirsadeghi. Were rnns all we needod? arXie preprint arXiv:2410.012 2024.
Leo Gao, Stella Biderman, Sid Black, Laurence Golding, Travis Hoppe, Charles Foster, Jason Phang, Horace He, Anish Thite, Noa Nabeshima, Shawn Presser, and Connor Leahy. The pile: An 800gb dataset of diverse text for language modeling, 2020.
William Gilpin. Chaos as an interpretable benchmark for forecasting and data-driven modelling, 2023a.
William Gilpin. Model scale versus domain knowledge in statistical forecasting of chaotic systems, 2023b.
William Gilpin. Generative learning for nonlinear dynamics. Nature Reviews Physics, 6(3):194-200 2024.
Xavier Gonzalez, Andrew Warrington, Jimmy T.H. Smith, and Scout W. Linderman. Towards Scalable and Stable Parallelization of Nonlinear RNNs. Advances in Neural Information Processing Systems, 37: 5817 5849, 2024.
Riccardo Grazzi, Julien Siems, Jörg KH Franke, Arber Zela, Frank Hutter, and Massimiliano Pontil. Un- locking slate-tracking in linear rans through negalive eigenvalues. arXiv preprinl arXiv:2411. (0or2/111.2599 2024.
Albert Gu and Tri Dao. Mamba: Linear-time sequence modeling with selective state spaces, arXio preprint arXi::231..00752 2023.
16

---

Published in Transactions on Machine Learning Research (09/2025)
Albert Gu, Karan Goel, and Christopher Ré. Efficiently modeling long sequences with structured state spaces. In International Conference on Learning Representations (TCLR), 2021.
Albert Gu, Karan Goel, Ankit Gupta, and Christopher Ré. On the parameterization and initialization of diagonal state space models. Advances in Neural Information Processing Systems, 35:35971 35983, 2022.
John L Gustafson and Isaac T Yonemoto. Beating floating point at its own game: Posit arithmetic. Super- computing frontiers and innovations, 4(2):71-86, 2017.
Sepp Hochreiter and Jürgen Schmidhuber. Long short-term memory. Neural computation, 9(8):1735-1780, 1997.
N.G. Kingsbury and P.J.W. Rayner. Digital filtering using logarithmic arithmetic. Electronics Letters, 7: 56 56-58, 1971. doi:
Toannis Kouretas and Vassilis Paliouras. Logarithmic number system for deep learning. In 2018 7th Inter- national Conference on Modern Circuits and Systems Technologies (MOCAST), pp. 1-4. IEEE, 2018.
Yann LeCun, D Touresky, Hinton, and T Sejnowski. A theoretical framework for back poockpopppggtio. In Proceedings of the 1988 connectionist models summer school, volume 1, pp. 21-28, 1988.
Yann LeCun, Corinna Cortes, and CJ Burges. Mnist handwritten digit database. ATT Labs [Online]. www.n..eemeemmmmmmmm coom/eddb/mann 2, 2010.
Available:
Yi Heng Lim, Qi Zhu, Joshua Selfridge, and Muhammad Firmansyah Kasim. Parallelizing non-linear se- quential models over the sequence length. In The Twelfth International Conference Learning Repre- sentations, 2024.
Winfried Lohmiller and Jean-Jacques E Slotine. On contraction analysis for non-linear systems. Automatica, 34(6):683-699 1998.
Jonas Mikhaeil, Zahra Monfared, and Daniel Durstewitz. On the difficulty of learning chaotic dynamics with rmns. Advances in neural information processing systems, 35:11297 11312, 2022.
Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio. On the difficulty of training recurrent neural networks. In International conference on machine learning, pp. 1310 1318. Pmlr, 2013.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Z. Yang, Zach DeVito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Sounith Chintala. Pytorch: An imperative style, higl--peffrmmnce deep learning library. CoRR, abs/1912.01703, 2019.
Arkady Pikovsky and Antonio Politi. Lyapunou exponents: tool to explore compler dynamics. Cambridge University Press, 2016.
Sheldon M Ross.
Stochastic processes. John Wiley & Sons, 1995.
Arnab Sanyal, Peter A Beerel, and Keith M Chugg. Neural network training with approximate logarithmic computations. In ICASSP 2020-2020 IEFE International Conference on Acoustics, Speech and Signal Processing (ICASSP), pp. 3122 3126. IEEE, 2020.
Jimmy T.H. Smith, Andrew Warrington, and Scott W. Linderman. Simplified state space layers for sequence modeling. In International Conference on Learning Representations (ICLR), 2023.
Anand Srinivasan and Jean-Jacques Slotine. Contracting differential equations in weighted banach spaces. Journal of Differential Equations, 344:203-229, 2023.
Steven H Strogatz. Nonlinear dynamics and chaos with student solutions monual: With applications to physics, biology, chemistry, and engineering. CRC press, 2018
17

---

Published in Transactions on Machine Learning Research (09/2025)
Philippe Tillet, H. T. Kung, and David Cox. Triton: intermediate language and compiler for tiled neural network compuiations. In Proceedings of the 3rd ACM SIGPLAN International Workshop on Machine
Learning and Programming Languages, MAPL 2019, pp. 10 10-19, 19, New York, NY, USA, 2019. Association 00145//311500832220
for Computing Machinery. ISBN 9781450367196. doi:
Songlin Yang, Bailin Wang, Yu Zhang, Yikang Shen, and Yoon Kim. Parallelizing lincar transformers with the delta rule over sequence length. In Proceedingg of NewrIPS, 2024.
18

---

Published in Transactions on Machine Learning Research (09/2025)
A Time to Estimate Spectrum of Lyapunov Exponents by System
multiple The plots below show time to estimate the spectrum of Lyapunov exponents sequentially, of time to estimate it in parallel with our algorithm, as we increase the number of time steps, for every dynamical system a dataset spanning multiple scientific disciplines (Gilpin 2023a bb. Each point is the mean seven runs on a single Nvidia GPU, In a few cases (most notably, "MacArthur)) ",, the improvement exccution time tapers off as we increase the number steps to and above 10 because parallel computation of all QR decompositions saturates the single GPU approximately 100% utilization.
Atrwve
www.eoiiieeeeee chyy
Amaadd
Se't Wettowcheeeeeeeo medi
Diodeytt
06% ryfadd
20"
Sohhen
DestefeerA..
DuPlyy
ExtendeCe
Frence
Pasfeerccee
Paredtrat Maggre
Focadhe
Spatioe
Spodi
Tyeetorr
Hyno
Hearlbing
Hpents
Oparies
I6C
Packtar
Spxits
Neath
Saatt
Seadte
Sealte
15
Necse pl
Tamrs
M6 Walhehraw
N00wag
107
10" 12"
157
19
307
107
10"
107
15'
127
107
10
10
12'
127
80
157 127 12"
19

---

Published in Transactions on Machine Learning Research (09/2025)
B
Derivation of Expression for Estimating Largest Lyapunov Exponent
If we express the sum of scalar logarithms in equation 23 as logarithm of scalar products, we obtain:
sell
LLE
log
u0
AtT
t=1
|ST||
1
|OT |uol
AIT
log
TT
L
x
|T
log
ArT
log || .rr....l J2. uoll
AfT
log
Σ (Jr )2
ArT
(29)
Norm computed
1
log
EJr ..... (JJ
2ATT
log
exp(log(Jr J2J(uO))²
2AIT
LSE (2 L.ME Jr. LMME(.... LMME(10g J2 LMMEE(og log (o)))))
24T7
LSE (2 LMMB(J,, LMME( LMME(J', LMMB(JI, w()))))
2AIT
LSE(2 SEN/PPPCCALLMMM)CCCM ... u6)).
24TT
20

---

Published in Transactions on Machine Learning Research (09/2025)
C Intuition behind Selective Resetting, with Step-by-Step Examples
Consider, for a moment, the following timo-variant non-diagonal lincar recurrence with biases,
X, = A,X,-1 + Br,
(30)
where X,, A,, B₁, and initial state Xo are all square matrices in Rdxd for (1.2.....n)) If we zero-out a particular transition matrix At, we "reset" the recurrence that step, with the corresponding Be becoming the new initial state. For example, given recurrence with 3 states,
1 = A1 + B1 // Xo is the initial state
2 = A2X + B2
(31)
X3 = A3X + B3,
if we set A2 - [0]axd, then B2 becomes the new initial state for subsequent states:
X1 X1 + Bı
X₂ ArXT B2 11 B. becomes the initial state for subsequent states
(32)
X3 = A₂X₂ + B3.
We can apply the same mechanism to reset linear recurrences without biases. Initialize all B to all-zeros, such that, initialization, we have the equivalent of a recurrence without biases:
X, = A.X, + BG // B,'s zeroed-out initialization
(33)
Now, at any step, we can reset this linear recurrence without biases, by zeroing-out that step's transition matrix and replacing the corresponding bias matrix with new initial stale of our choice. For example, a sequence n = 3 states with all B''s initialized to all-zeros, is:
X1 A1Xo + Bi
Axoo
1/ No is the initial state
X2 A2X1 + W2
A2X1
(34)
X3 = AzX₂ + DJ/ = =A3X2.
If we set, say, B2 - R(A2X1), then A2 - (0)axd, where R. Raddd b Rdxd is function of our choice,
obtain the following sequence of states:
X1 = A1 + B = Ai
X₂ = ArXI + B2 = B2
1/ Bg becomes the new initial state for all subsequent states
(35)
X3 = A&X2 + Bi = A3X2.
Now, consider parallel prefix scan, applied right-to-left, over the unmodified recurrence.
21

---

Published in Transactions on Machine Learning Research (09/2025)
C.1
Parallel Prefix Scan without Selective Resets
When we apply parallel prelix sean, right-to-left, over the ummodified recurrence, we have:
Parallel Prefix Scan:
A3
A2
A
Xo
Input states:
bias states initialized all-zeros
B3
B2
B1
Bo
A3
A2
Ao
Xo
Parallel step :
B2 + B3
update
B3
B + B
update
Bo
Ay (
A Xp)
A1 A Bo + B1
Xo Bu
Parallel step 2:
A, + B1) + Ag + B3
A2(A, B1) + B2
update
A1A2A1Xo {0]axd
AoA1Xo (0)dxd
A [0)dxd
Xo {0]dxd
Output states:
// Bj all-zeros
Sums:
Xo
Ao
X
X
|/ same as recurrence without blases
(36)
That is, we obtain X; = A,X1-1, equal to the original recurrence without biases.
C.2 Parallel Prefix Scan with One Selective Reset
If we selectively reset, say; the second interim compound state, before parallel step 2, we have:
Parallel Prefix Scan, Selectively Reset before Parallel Step
Aa
A2
A
X
Input states :
// bias states initialized to all-zeros
B3
B2
B
Bo
Aq.A2
A2
A Xo A Bu + B
Xo Bu
Parallel step 1:
Ag B2 + B3
update
B2
update
Ao A Bo B1
{0)dxd
Selective reset:
Replace
with
(37)
R(A,Xo)
Ag Ag d0]dxd) A3.A2 R(A,Xo) + A3 + BB
A2(10)dxd) Ag Xo) + BB
{0]dxd
Xn
Parallel step 2:
update
R(A,Xo)
Bo
{0]axd AsA2 R(AIXo)
{0)axd
{o)axa R(A1Xo)
Xo {O]axd
Output states:
11 Be's all-zeros
A2 R(A
Sums
AsA2 R(A,X()
A R(A
R(A,
X
// modified linear recurrence
That is, we obtain a modified sequence in which step 2's state was selectively reset, changing its value from X, to Xp), compounded in all subsequent states, as we computed all states in parallel.
22

---

Published in Transactions on Machine Learning Research (09/2025)
D Quantitative Comparisons of Our Implementation to Floating-Poine Point Formats
We compare Complex128 and Complex64 GOOMs, respectively, to Float64 and Float32, the two float formats with greatest precision and dynamic range currently supported on Nvidia GPUs. The comparisons, shown in the pages that follow, are valid only for our implementation, not for GOOMs in general.
Magnitude of Errors We compare the maguilude of errors for one- and two-argument scalar functions that can be used to compose many other functions, and also for a representative matrix product. The one- argument scalar functions are: reciprocal y = 1/z1, square root y = vz, square y= = natural logarithm y logæ, x, and exponential y e". The two-argument scalar functions are: addition = x + y (and implicitly, subtraction, since y = + -y) and scalar product = zy (and implicitly, division, since z/v = ry-1). The matrix product is of two 1024 x 1024 matrices with elements independently sampled from (N 1), representative of typical matrix products in deep learning models. We apply all one- and two-argument functions except the exponential function over a range of values spanning the approximate number of decimal digits to which each float format is precise: 10-15 to 1015 for Float64, and 10 to 10 for Float32. We apply the exponential function over a range of values spanning 10-5 to 10. For one-argument functions, the range consists of 1M input values, equally spaced in the domain of decimal digits. For two- argument functions, each argument's range consists of 10,000 input values, equally spaced in the domain of decimal digits, such that we apply cach two-argument function to 100M different pairs of input values. For all one- and (wo-argument functions, we measure the magnitude of absolute errors in base 10 (i.e., the number of decimal digits of error) in relation to the same operation executed over Float128, a float format with greater procision and dynamic range than Float64 and Float32. For the matrix product, we measure the error normalized by the Frobenius norm of the matrix product in relation to the same operation executed over Float128. We execute all operations on a recent Nvidia GPU, except for operations and comparisons over F28at128, because Nvidia GPUs currently do not support it. We implement all transformations over complex-typed GOOMs so they accept floats as inputs, internally map them to complex-typed GOOMs, execute all computations over complex-typed GOOMs, map them back to floats, and, finally, return those floats as outputs. Mapping complex-typed GOOMs to floats impacts precision but is necessary for comparison to Float128. We find that the magnitude of errors varies from approximately the same to only slightly more, notwithstanding the impact of mapping complex-tye typed GOOMs to floats.
Running Time We compare rumning time for every one- and two-argument scalar function, as we apply it to batch of 100M input samples, each independently drawn from (4(0,1) 1), processing each batch in parallel on a recent Nvidia GPU. We also compare running time for the representative matrix product. We apply each transformation 30 times, each time to a different batch, and compute the mean running time over all batches. We report the mean running time of each transformation over Complex and Complex64 GOOMs as a multiple of the mean running time for Float64 and Float32, respectively. We implement every transformation over complex-typed GOOMs to accept complex-typer GOOMs as inputs and return complex-yppd GOOMs as outputs, regardless of whether doing so penalizes running time. For example, we implement scalar addition over complex-typed GOOMs as log(e" e²', not as ex + e² significantly impacting running time due to the creation of float tensors for storing interim exponentiated values and complex tensors for storing output values. On the other hand, our implementation of natural logarithm incurs zero running time, because complex-typed GOOMs are already natural logarithms. We find that for most operations, the running time of our implementation is approximately twice that of floats.
We compare memory use for every one- and two-argument. scalar function, as we apply it to batch of 100M input samples, each independently drawn from 2((0, 1), processing each batch in parallel recent Nvidia GPU. We also compare memory use for the representative matrix product. For both complex-typed GOOMs and floats, measure peak memory allocated, including creation of input, interim, and output tensors, as well as PyTorch overhead, if any. We report peak memory allocated over Complex122 and Complex64 GOOMs as a multiple of peak memory allocated over Float64 and Float32, respectively. We find that peak memory allocated is typically at least twice that of floats, but sometimes it can be less.
Memory Use
For all comparisons, see the pages that follow. To replicate the comparisons, see our published source code.
23

---

Published in Transactions on Machine Learning Research (09/2025)
D.1 Complex128 GOOMs versus Float64
We show the comparisons of Complex128 GOOMs to Float64 first, because they likely important to researchers and practitioners interested in quantitative comparisons of precision.
D.1.1
Magnitude of Errors on Reciprocals
Magnitude of Error versus Float128 for Reciprocals, y = 1/x Complex128 GOOM
Float64
-5
-15
l0g1o ly
-20
-25
-35
-15
-10
-5
0
5
10
15
-15
-10
-5
0
5
10
15
y:
log1oX
log10
Tested Value
D.1.2 Magnitude of Errors Square Roots
Magnitude Error versus Float128 for Square Roots, y = =
Complex128 GOOM
Float64
-7.5
110.0
-12.5
-15.0
log10 ly
-17.5
-20.0
-22.5
-25.0
-5
-15
-10
15
Plat128
log1oX
logio
Tested Value
24

---

Published in Transactions on Machine Learning Research (09/2025)
D.1.3
Magnitude of Errors Squares
Magnitude of Error versus Float128 for Squares, y = x² Complex128 GOOM
Float64
10
log ly
-30
-15 -10
-5
5
10
15
-15
-10
-5
5
10
15
y: Flas1288
log10 X
log10 X
y: Tested Value
D.1.4
Magnitude of Errors Natural Logarithms
Magnitude of Error versus Float128 for Natural Logarithms, y log x
Complex128 GOOM
Float64
-15
-16
logioly ly yl
-17
-18
-19
-15
-10
-5
10
15
-15
-10
-5
5
10
15
Plat128
log10
log10
Tested Value
25

---

Published in Transactions on Machine Learning Research (09/2025)
D.1.5
Magnitude of Errors Exponentials
Magnitude of Error versus Float128 for Exponentials, y = Complex128 GOOM
Float64
-10
-12
10g:0 by yl
-14
-18
-
-3
-2
-1
0
1
-5
-4
-3
-2
-1
0
y: Ros1128
log1o
log10 x
2: Texced Value
D.1.6
Magnitude of Errors on Scalar Addition/SSubraationn
Magnitude of Error versus Float128 for Scalar Addition, z x +
Complex128 GOOM
Float64
15
-5
5
-10
logioy
0
log:9
-20
-5
-10
330
-15
-15
-10
-5
10
15
-15
-10
-5
0
5
10
15
1
1og X
Fleat:20
26

---

Published in Transactions on Machine Learning Research (09/2025)
D.1.7
Magnitude of Errors on Scalar Product /Diiisinn
Magnitude of Error versus Float128 for Scalar Product, = xy
Complex128 GOOM
Float64
15
10
5
log1o
-10
0
10g
20
-5
-3
40
-15
-15
-10
-5
-15
-10
-5
5
10
1og1oX
logro
Fbatt 28
Tested Va ue
D.1.8
Normalized Errors on Representative Matrix Product
Histogram of Normalized Errors versus Float128 Matrix Product, - where X, are 1024x1024 Matrices with Elements Sampled from ./ 1)
Complex128 GOOM
Float64
10
10-18
-5 x 10-10
-1
2500
5000
7500
10000
12500
15000
2500
5000
7500
10000
12500
15000
2: Deat125
Number of Elements
Number of Elements
2: Tested Value
27

---

Published in Transactions on Machine Learning Research (09/2025)
D.1.9
Execution Times
Execution Time on an Nvidia GPU, Complex128 GOOM as Multiple of Float64
Reciprocals
Square Roots
Squares
Natural Logarithms
Exponentials
Scalar Addition
Scalar Product
Matrix Product
10
Execution Time, Complex128 GOOM as Multiple of Float64 (Mean of 30 Runs)
100M
Matrix product: 1004 mats
Mean execution times are reported by torch Norch.utila wtw.lls..bnchaa.... Timer. To the extent possible, apply transformations in-place, to minimize the impact of new memory allocations execution time. Complex128 GOOMs are already natural logarithms, so obtaining them incurs no computation.
D.1.10
Peak Memory Allocated
Peak Memory Allocated on an Nvidia GPU, Complex128 GOOM a Multiple Float64
Reciprocals
Square Roots
Squares
Natural Logarithms
Exponentials
Scalar Addition
Scalar Product
Matrix Product
0.0
0.5
1.0
1.5
2.0
2.5
3.0
Peak Memory Allocated, Complex128 GOOM Multiple of Float64 (Including Memory Allocated to Input, Interim, and Output Tensors)
Scalor fu nctions: 1030 in peral
Memory allocation figures are as reported by torch. cuda salaaaammory..... max memory allocated. the extent possible, do not apply transformations in-place, to allocate memory for input, interim, and output tensors. Memory by the matrix product Complex123 GOOMs is for initial implementation of LMME.
28

---

Published in Transactions on Machine Learning Research (09/2025)
D.2 Complex64 GOOMs versus Float32
We show the comparisons of Complex64 GOOMs to Float32 in the same order the comparisons of Com- plex128 GOOMs to Float64, for easier cross-reference.
D.2.1
Magnitude of Errors on Reciprocals
Magnitude of Error versus Float128 for Reciprocals, y = 1/x
Complex64 GOOM
Float32
0.0
-25
-50
-7.5
log10
-10.0
-12.5
-15.0
-17.5
-2
0
6
-4
-2
6
V:
log1oX
lagio
Tested Value
D.2.2
Magnitude of Errors Square Roots
Magnitude Error versus Float128 for Square Roots, y = =
Complex64 GOOM
Float32
-4
-8
log10 ly y/
-12
-14
-2
-4
-2
Poat128 Tested Value
log10 X
log10
29

---

Published in Transactions on Machine Learning Research (09/2025)
D.2.3
Magnitude of Errors on Squares
Magnitude of Error versus Float128 for Squares, y = x2 Complex64 GOOM
Float32
-5
lagoo ly -y
-10
-15
-20
-25
4
-2
2
4
6
-4
-2
2
y: Flas1288
logi0 X
log10 X
V: Tested Value
D.2.4
Magnitude of Errors Natural Logarithms
Magnitude of Error versus Float128 for Natural Logarithms, y log x
Complex64 GOOM
Float32
-6
-7
-8
loggaly ly yl
-9
-10
-12
-2
-2
2
Poat128 Tested Value
log1o X
log1o
30

---

Published in Transactions on Machine Learning Research (09/2025)
D.2.5
Magnitude of Errors Exponentials
Magnitude of Error versus Float128 for Exponentials, y == Complex64 GOOM
Float32
-2
-4
-6
log:o by -y
-10
-12
14
-5
-3
-2
-1
0
1
-5
-4
-3
-2
0
7: Ros1128
lag10 x
log10 x
2: Texced Value
D.2.6
Magnitude of Errors on Scalar Addition/SSubraationn
Magnitude of Error versus Float128 for Scalar Addition, z +y
Complex64 GOOM
Float32
0.0
-2.5
-5.0
logsay
-7.5
-10.0
10g10
-2
-12.5
-15.0
4
-17.5
-6
6
-4
-2
0
6
6
-2
0
2
1
log1oX
Flbat::0
31

---

Published in Transactions on Machine Learning Research (09/2025)
D.2.7
Magnitude of Errors Scalar Product /Diiisinn
Magnitude of Error Float128 for Scalar Product, xy
Complex64 GOOM
Float332
5
l0g:oy
-5
109
-10
-4
-6
-6
-4
-2
0
2
6
6
-4
-2
2
6
10
logio
Fibat 20
Tested Va
D.2.8
Normalized Errors on Representative Matrix Product
Histogram of Normalized Errors versus Float128 Matrix Product, Z= where X, are 1024x1024 Matrices with Elements Sampled from ./ 1)
Complex64 GOOM
Float32
10
100
-5 x 10-10
-10-9
2000
6000
8000
10000
2000
4000
6000
8000
10000
2: Deat125
Number of Elements
Number of Elements
2: Tested Valuo
32

---

Published in Transactions on Machine Learning Research (09/2025)
D.2.9
Execution Times
Execution Time on an Nvidia GPU, Complex64 GOOM as Multiple of Float32
Reciprocals
Square Roots
Squares
Natural Logarithms
Exponentials
Scalar Addition
Scalar Product
Matrix Product
2
Execution Time, Complex64 GOOM as Multiple of Float32 (Mean of 30 Runs)
100M
Matrix product: 1004 mats
Mean execution times are reported by torch.utin wut 1s. benchmaa...............mm Timer. To the extent possible, apply transformations in-place, to minimize the impact of new memory allocations on execution time. Complex64 GOOMs are already natural logarithms, so obtaining them incurs no computation.
D.2.10
Peak Memory Allocated
Peak Memory Allocated an Nvidia GPU, Complex64 GOOM Multiple of Float32
Reciprocals
Square Roots
Squares
Natural Logarithms
Exponentials
Scalar Addition
Scalar Product
Matrix Product
0.0
0.5
1.0
1.5
2.0
2.5
3.0
3.5
Peak Memory Allocated, Complex64 GOOM as Multiple Float32 (Including Memory Allocated to Input, Interim, and Output Tensors)
Scalor functions: 10001 in panal
Memory allocation figures are as reported by torch. monda.memory.m.. max memory allocated. To the extent possible, do not apply transformations in-place, to allocate memory for input, interim, and output tensors. Memory by the matrix product Complex64 GOOMs is for initial implementation of LMME.
33


Reviews:
[{"review_id": "rK5qBO4NbA", "reviewer": "Reviewer_Tz4F", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n1. Novel idea: Encoding real numbers as complex logarithms seems to be impactful.\n\n2. Breadth of applications: Demonstrated in numerical linear algebra, dynamical systems, and deep learning.\n\n3. Practical implementation: GOOMs are implemented in PyTorch with GPU support and autograd compatibility.\n\n4. Parallelization aspect: Combining GOOMs with prefix scans and selective resetting is technically interesting.\n\nWeaknesses:\n\n1. Positioning vs. prior work: The distinction from logarithmic number systems (LNS) and prior log-domain methods is not sufficiently clear.\n\n2. Empirical validation: Experiments are mostly small-scale or synthetic; comparisons against established stabilization techniques (e.g., normalization, gradient clipping, high-precision arithmetic) are missing.\n\n3. Performance: The compromise implementation of log-matrix-multiplication-exp (LMME) incurs overhead (~2× slower) without thorough benchmarking.\n\n4. Clarity: Parts of the exposition (e.g., selective resetting, definition of GOOMs) are abstract and mathematically heavy, which may limit accessibility.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper presents Generalized Orders of Magnitude (GOOMs), a complex-logarithmic representation of real numbers designed to enable stable computation across dynamic ranges beyond floating-point limits. The contributions are: (1) formalization of GOOMs and demonstration that floating-point numbers are a special case; (2) implementation of GOOMs in PyTorch with GPU/autograd support; (3) introduction of a selective-resetting method for parallel prefix scans; and (4) demonstration of applications in long chains of matrix products, Lyapunov exponent estimation, and recurrent neural networks.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. Clarify Relation to Prior Work: The paper could more clearly articulate how GOOMs differ from and extend traditional logarithmic number systems (LNS) and log-domain arithmetic methods. At present, the novelty may appear incremental without explicit comparisons and positioning.\n\n2. Expand Empirical Validation: Current experiments are largely synthetic or small-scale. It would be important to include evaluations on real-world machine learning or scientific computing tasks, and to compare against standard stabilization techniques (e.g., gradient clipping, normalization methods in RNNs, extended precision or log-space arithmetic).\n\n3. Performance Benchmarks: Provide more detailed benchmarks quantifying the computational and memory overheads of GOOMs relative to conventional floating-point methods. For example: FLOPs, runtime, GPU utilization, and scalability beyond the toy settings presented.\n\n4. Improve Exposition and Accessibility: Sections on the definition of GOOMs and the selective-resetting method could benefit from clearer intuition, figures, or worked examples earlier in the main text. The appendix is helpful, but the core paper would be more readable if some of that material were integrated.\n\n5. Applications Beyond Demonstrations: If possible, demonstrate GOOMs on a larger-scale deep learning task (e.g., a benchmark language model or transformer setting), or a scientific computing workflow, to highlight practical impact and not just conceptual novelty.", "broader_impact_concerns": "This work is primarily technical in nature, focusing on numerical representations and stability in high-dynamic-range computation. As such, it does not raise immediate ethical concerns typical of machine learning applications (e.g., bias, fairness, misuse of generative models).\n\nThe main broader-impact considerations relate to:\n\n1. Computational Efficiency and Resource Use: If GOOMs incur significant overhead compared to floating-point arithmetic, widespread adoption could lead to higher energy consumption and reduced efficiency in large-scale ML training. The paper may need to briefly discuss these trade-offs.\n\n2. Accessibility and Usability: Introducing a new numerical representation could affect software and hardware ecosystems. While this work remains at the software-level, readers would benefit from a reflection on whether adoption poses barriers for practitioners unfamiliar with complex arithmetic.\n\nOverall, there are no major ethical risks, but a short broader-impact statement acknowledging the potential implications for compute efficiency, sustainability, and adoption costs would make the submission more complete.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "MoVLFi8NtJ", "reviewer": "Reviewer_o1nF", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strengths\n\n- Clear presentation\n- Tackles a significant problem of very deep or recurrent architectures trained over long time horizons\n- Introduces a parallelized algorithm for computing Lyapunov spectra, based on the author's resetting method\n- Implementation in pytorch\n\n## Weaknesses\n\n- Almost no review of previous related literature\n- Numerical experiments could be strengthened by demonstrating small numerical error, beyond avoidance of under/overflow\n- Clearer delineation of limitations would make the method easier to re-use\n\nOverall, I believe the manuscript makes a significant contribution, is clearly presented, and I did not identify technical issues (though see restriction below). The main weakness that should be addressed is the lack of embedding in the state of previous related literature, since the use of the log domain to avoid float over/underflow is a well-known approach. \n\nPlease note that I am not an expert in low-level numerical computing, and have therefore reviewed the manuscript on a less detailed level.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors propose a software-level method representing real numbers as complex logarithms (the subset of the complex plane that exponentiates to ℝ), which is inspired from earlier logarithmic number systems. They demonstrate how this allows for numerically stable computation across a larger dynamic range than float32/64, with applications including long matrix product chains such as in Lyapunov spectrum estimation and and training RNNs. The implementation is integrated into PyTorch and exploits GPU parallelism via prefix scans. In addition, the authors provide a novel parallelized algorithm for computing the Lyapunov spectrum of high-dimensional systems.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- The authors mention that their approach builds on prior work on logarithmic number systems (LNS), going back to the work by Kingsbury and Rayner (1971). Also the log-sum-exp trick has been standard in this literature. As a non-expert reader, my impression is that the conceptual innovation here to use the multiple of $\\pi i$ to represent the sign, however the implications of this change compared to manually representing the sign where not fully obvious to me. Please clarify, and in particular provide a more detailed review of previous and related work on LNS or similar approaches.\n- To a reader not familiar with the state of the art in estimating Lyapunov spectra, the introduction and discussion imply that previous approaches where always sequential instead of parallel as per the proposed method. In section 4.2 one reference is given regarding this. Could the authors add a short paragraph giving an overview of previous methods?\n- Section 2: What kind of transformations leveraging the structure of the complex plane (e.g., rotations in $\\mathbb{C}$) are allowed that would not break exponentiation to $\\mathbb{R}$?\n- The experiments demonstrate that the computation does not run into over- or underflow, but there is no evidence of the precision of the results. Is it feasible and interesting to add measurements of numerical error for a toy-task example by comparing to analytically computed values? For example, this could be a simple toy task of matrix exponentiation, or comparison of the computed Lyapunov spectrum to an analytically known spectrum.\n- The manuscript clearly explains why the logarithmic numbers have a higher dynamic range, and that the representable numbers when using floats as the GOOM exponents are gracefully distributed across the higher range. It is also clear that multiplication will be numerically stable with GOOMs for (very) large numbers, and that the experiments demonstrate the practical success on the example applications (although see the question on precision above). However, it appears that using the system must also come with costs, which are currently not clearly articulated. Firstly, the difficulty of performing addition and in how far it is overcome by the log-sum-exp trick could be explained more clearly, and existing literature studying this trade-off should be mentioned. Secondly, are there operations where the numerical error due to floating point precision will be greater for float+GOOM than for floats themselves? As a non-expert reader I came away from the manuscript with the impression that GOOMs could be a good tool, but that I am not sure where this tool may fail, which is of course detrimental to adoption. Giving clear guidelines or references on such limitations would therefore improve the impact of the paper. I leave the extent to which these questions should be answered by additional analysis or improved discussion to the authors.", "broader_impact_concerns": "/", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "8VOoIX9LSv", "reviewer": "Reviewer_djDE", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength:\n* The proposed method is easy to implement, intuitive, and effective in expanding the dynamic range of traditional floating-point numbers.\n* The authors conduct various experiments to test the effectiveness of their method. They also design and implement an efficient parallel algorithm for computing the Lyapunov spectra.\n\nWeakness:\n* The idea of using log to handle large number multiplications is not new (e.g., [1]).\n* The proposed method is not very efficient, according to my understanding. While it simplifies multiplication, it makes addition much more complicated, considering that taking the exp/log is relatively slow. Additionally, half of the storage (the imaginary part) is used solely to store a binary message (the sign), resulting in a significant waste of memory.\n* No running time/memory usage/stability analysis results are reported. It is not clear how the proposed method compares to floating-point numbers in these aspects.\n\n[1] Parhami B. Computing with logarithmic number system arithmetic: Implementation methods and performance benefits. Computers & Electrical Engineering. 2020 Oct 1;87:106800.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces GOOM, an alternative for storing real numbers as complex-valued representations. The real part is used to store the log-magnitude, and the imaginary part is used to store the sign. GOOM extends the dynamic range of traditional floating-point number formats like float32/64. Using GOOM, the authors implement the computation of matrix multiplications, Lyapunov spectra, and neural network training.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "* Add further experiments on running time/memory usage/relative error comparisons between GOOM and traditional floating-point numbers.\n* Show how GOOM is compared to alternative extended-range formats, such as Posit Arithmetic[2].\n\n[2] Gustafson JL, Yonemoto IT. Beating floating point at its own game: Posit arithmetic. Supercomputing frontiers and innovations. 2017 Apr 25;4(2):71-86.", "broader_impact_concerns": "No concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)