Pith. sign in

REVIEW 4 major objections 5 minor 16 references

A Systolic Array Architecture for Nonlinear Activation Functions and Softmax Computation using Chebyshev Polynomials

T0 review · 4 major / 5 minor · reviewed 2026-08-06 · deepseek-v4-flash

Pith's one-line read This paper claims that one fixed-point systolic array can compute tanh, sigmoid, ReLU, and softmax by evaluating Chebyshev polynomial approximations with Horner's method, achieving up to 71% lower tanh error than a CORDIC baseline while…

desk verdict A clean Horner-based systolic activation unit with a genuine resource-sharing idea, but the headline softmax accuracy gains rest on an unstated KL direction and an assumed baseline; the engineering numbers are solid. read the letter →

arxiv 2608.04734 v1 pith:OFU7C6C2 submitted 2026-08-05 cs.AR

classification cs.AR
keywords ArrayandvectorprocessorsChebyshevapproximationtheoryNeuralnetsSpecial-purposehardwaresoftmaxaccelerationactivationfunctionsfixed-pointarithmeticHorner'smethod
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

This paper claims that a systolic array—a grid of multiply-add cells normally used for matrix multiplication in neural accelerators—can double as a universal nonlinearity unit. By evaluating Chebyshev polynomial approximations with Horner's method, the same hardware computes tanh, sigmoid, and ReLU for individual neurons and softmax for an entire output vector. The paper reports that this shared unit cuts tanh mean absolute error by up to 71% versus a feature-matched CORDIC design, uses 4.6% less area and 5.1% less power, and reduces softmax KL divergence by 44.6% versus CORDIC and 79.0% versus a piecewise-linear baseline at input dimension 256. If these numbers hold, accelerators for transformers and classifiers, where softmax is a significant runtime share, could fold activation and softmax into one compact fixed-point datapath.

What carries the argument

The load-bearing mechanism is Horner's method for polynomial evaluation combined with Chebyshev polynomial approximation. Horner's method rewrites a degree-n polynomial as a0 + x(a1 + x(a2 + ... x(an))), a chain of multiply-add steps. A systolic array's processing elements can each perform one multiply-add per cycle, so a column of n+1 PEs evaluates a degree-n polynomial at one input, with coefficients a_i propagating through the array and x held constant down the column. Chebyshev polynomials supply near-minimax approximations, so low-degree polynomials (degree 9 for tanh, degree 5 for sigmoid, degree 3 for e^x) suffice for the target accuracy. For softmax, the shifted-softmax identity confines the exponential inputs to (-inf,0], making one e^x polynomial valid, and a hardware divider computes the reciprocal of the sum; clipping of e^x below x_min removes negligible probabilities while preserving the ranking and ratios of the survivors.

What would settle it

Replace the assumed ONE-SA softmax baseline with an implementation that uses a hardware divider and a more accurate e^x approximation, and re-measure the maximum KL divergence at N=256; if the gap narrows below the reported 79% or ONE-SA's KL drops below 3.252, the central softmax advantage would not hold. Independently, running the CORDIC baseline with the range-extension method referenced in [16] and re-measuring tanh MAE on [-4,4] would test whether the 71% error reduction persists.

Watch

Extended reading notes

Core claim

The central discovery is that a polynomial approximation evaluated on a systolic array can serve as a universal nonlinearity engine. Chebyshev polynomials of the first kind, computed in power-series form via Horner's method, map naturally onto the array's per-PE multiply-add: each column evaluates the polynomial for one input, with coefficients streaming down the array. The same hardware computes tanh and sigmoid directly, ReLU with a comparator, and softmax by approximating e^x on shifted inputs, summing the exponentials, and dividing. A quantization-aware, two-step search over polynomial degree, bitwidth, and clipping range selects a unified Q3.12 fixed-point format that balances the needs of all three functions. The measured result is a 10-row, 8-column systolic array that beats a feature-matched 8-unit CORDIC design in mean and max absolute error for tanh and sigmoid, reduces tanh/sigmoid latency by 16.2% and softmax latency by 14.3%, and uses less area and power for the same 8-wide parallelism.

Load-bearing premise

The headline error reductions rest on comparing against a CORDIC baseline that is allowed to diverge outside its convergence range and against an assumed piecewise-linear softmax baseline for ONE-SA, so the measured gaps partly reflect those baseline choices rather than the polynomial method itself.

Editorial extensions

If this is right

  • A single fixed-point systolic array can replace separate activation and softmax units in a neural accelerator, cutting activation-area by 4.6% and power by 5.1% compared with a CORDIC-based unit at equal 8-wide parallelism.
  • tanh and sigmoid can be computed with up to 71% and 41% lower mean absolute error than a feature-matched CORDIC design, with latency per output reduced by 16.2% for both functions.
  • Softmax KL divergence drops by 44.6% versus CORDIC and by 79.0% versus a piecewise-linear baseline at input dimension N=256, with row-sum error 0.021 versus 1.021 for the piecewise-linear baseline.
  • The same architecture extends to other univariate functions by changing hard-coded polynomial coefficients, since the array structure and fixed-point datapath are function-agnostic.
  • Because the polynomial degree and the number of array columns can be chosen independently, the design can trade approximation accuracy against parallelism and area per application.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • Beyond the paper's stated results, the error gap over CORDIC would likely narrow if the baseline included the range-extension method cited in the paper, so the 71% figure is best read as against an unextended CORDIC unit rather than all CORDIC designs.
  • The 79.0% softmax improvement over ONE-SA depends on the paper's stated assumption that ONE-SA computes softmax with a piecewise-linear reciprocal; if ONE-SA's actual softmax is more accurate, the KL ranking at N=256 could change.
  • The same array could be applied to other elementwise nonlinearities in modern networks (e.g., GELU, SiLU) or to normalization statistics, because the array only needs new coefficient sets and clipping ranges; the paper does not explore these extensions.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

4 major / 5 minor

Summary. The paper proposes a systolic-array-based activation unit that evaluates Chebyshev-polynomial approximations of tanh, sigmoid, and e^x using Horner's method, with a comparator for ReLU and a four-step shifted-softmax datapath. Coefficients are selected by a quantization-aware grid search over polynomial degree, bitwidth, and clipping range. The authors compare their implementation against a CORDIC baseline and an assumed ONE-SA-style piecewise-linear baseline, reporting lower mean absolute error for tanh/sigmoid, lower KL divergence for softmax at N=256, and lower area/power for the 8-register variant.

Significance. The architecture itself is sensible and has genuine value: it reuses the systolic-array datapath for Horner evaluation, supports multiple univariate functions and softmax with one multiplier/adder chain, and the two-step quantization-aware coefficient search is a systematic way to balance degree, bitwidth, clipping range, and hardware cost. The synthesis results are presented for a standard PDK with clear area, power, and latency tables. If the baseline and metric issues identified below are repaired, the paper would be a useful contribution to low-cost activation and softmax units. At present, however, the headline numerical claims overstate what is demonstrated: the tanh/sigmoid error reduction depends on a CORDIC baseline without range expansion, the softmax KL reduction depends on an assumed ONE-SA implementation and on an unstated divergence direction that rewards clipping, and the area/power savings apply to a configuration that cannot compute the N=256 softmax results.

major comments (4)
  1. [§IV-B, Table II] The reported finite KL values require clarification of the divergence direction. Since Section III-C clips the exponential approximation to zero below x_min = -5, the approximate softmax P has zero probability for any class whose shifted input is below -5; for the test distributions N(0, 4) and U(-6, 6) with N=256 such classes occur with probability one. For the true softmax Q, Q_i > 0 for all classes, so D_KL(Q || P) is infinite. The finite entries in Table II therefore imply that the authors computed D_KL(P || Q), but this direction is not stated in Section IV-B, and it assigns zero contribution to exactly the classes that the approximation drops. The comparison with CORDIC and ONE-SA is asymmetric because those baselines do not perform equivalent support truncation. Please report the direction explicitly, justify it, provide the reverse direction or a smoothed alternative (for example, D_KL with an epsilon floor), and report the number of clipped classes. Because Section III-D uses KL divergence to select the e^x polynomial, the choice of degree 3 with x_min = -5 may itself be an artifact of this direction.
  2. [§IV-B, Fig. 5 and Table I] The CORDIC baseline is used outside its native convergence range without range expansion. The text states that the CORDIC approximation stops converging around |x| > 1.1182 for tanh and outside |x| > 2.236 for sigmoid, yet mean and maximum AE are reported over [-4, 4]. The claimed 71% and 41% MAE reductions are therefore partly an artifact of the baseline diverging in regions that a practical CORDIC activation unit would handle by range expansion or by a separate saturation path. Either include a range-expanded CORDIC baseline (for example, using the method of [16]) or restrict the comparison to the CORDIC convergence interval and report the saturation/clipping behavior separately.
  3. [§IV-B, Table II and Abstract] The 79.0% improvement over ONE-SA is based on an assumed implementation, not on ONE-SA as published: the paper states, 'we assume a calculation method similar to our 4-step approach from Sec. III-C, but using a piecewise linear approximation of 1/x instead of a hardware divider.' This is a transparent caveat, but it makes the abstract's claim of lower KL divergence 'compared to ... a piecewise-linear approximation' and the conclusion's statement about outperforming both CORDIC and ONE-SA too strong. Please rephrase the claim as a comparison against an assumed ONE-SA-style baseline, and if possible evaluate the actual ONE-SA softmax or provide a sensitivity analysis over plausible ONE-SA arithmetic choices.
  4. [Abstract and §IV-C, Table III] The 4.6% area and 5.1% power savings are for the 8-parallel, 8-register design, which cannot compute the N=256 softmax results in Table II because the maximum softmax vector length is limited by the number of intermediate registers. The 256-register designs that support N=256 have higher area and power than the CORDIC baseline (for example, 106192 um2 and 0.9707 mW versus 89719 um2 and 0.8646 mW for the 8-column 256-register variant). The paper should state explicitly that no single configuration simultaneously realizes the area/power savings and the N=256 softmax error, or it should present a configuration that does.
minor comments (5)
  1. [§IV-A] Please report the final quantized polynomial coefficients for the selected tanh, sigmoid, and e^x approximations; the paper gives only the illustrative low-degree examples in Eqs. (1)-(2), which are not necessarily the final coefficients used in the evaluation.
  2. [§III-C] The sentence 'For these non-clipped inputs, both their ranking and the ratios between their probabilities are preserved' should specify that this holds for the approximate exponentials after clipping, not relative to the exact softmax distribution.
  3. [Fig. 5] The y-axis spans many orders of magnitude and the Chebyshev curve is very close to zero in parts of the plotted range; adding markers or a zoomed inset would help readers verify the reported mean and maximum AE values.
  4. [§IV-B, Table II] The label 'ONE-SA (0.25)' could be confused with the actual ONE-SA design; consider renaming it to 'assumed ONE-SA variant with segment width 0.25' throughout.
  5. [§IV-B] Please define Row-Sum Error precisely (for example, max over samples of |sum_i p_i - 1|) and state the total number of samples over which the maximum is taken beyond the phrase '128 samples per distribution'.

Circularity Check

0 steps flagged · score 2.0 of 10

No significant circularity: Chebyshev coefficients are fitted approximations and the reported errors are evaluations, not re-statements; the self-cited CORDIC baseline and unstated KL-divergence direction are caveats, not circular reductions.

full rationale

The claimed derivation is not circular. The Chebyshev coefficients in (1)-(2) are obtained by interpolating tanh and e^x at Chebyshev nodes (Sec. II), and the accuracy numbers in Tables I and II are computed evaluations of those polynomials against the target functions/softmax, not re-statements of the fitted coefficients. The architecture is a standard Horner-method systolic dataflow; nothing in Sec. III defines an output in terms of the target metric. Two caveats explain the score rather than a circular step. First, the CORDIC baseline is taken from Ref. [9], which shares a co-author with this paper and is extended in-house; this is a benchmark-choice risk (the baseline stops converging near |x|>1.1182 for tanh), but it is an implemented, synthesizable design, not an unverified self-citation invoked to forbid alternatives. Second, Sec. IV-B reports finite KL values without stating the direction; since Sec. III-C clips e^x to zero below x_min=-5, D_KL(true||approx) would be infinite for some test distributions, so the finite table implies D_KL(approx||true). That asymmetry can inflate the claimed softmax improvement over baselines that do not truncate support, but it is an evaluation-metric concern, not a reduction of the prediction to its inputs. The assumed ONE-SA softmax (PLA-based 1/x) is explicitly labeled an assumption, so it does not masquerade as an external result. Verdict: no significant circularity.

Assumptions & free parameters 4 free parameters · 5 assumptions · 0 invented entities

The central claims rest on fitted polynomial coefficients and clipping bounds chosen by search, on the assumption that quantization-aware simulation matches the synthesized 16-bit hardware, and on two baseline modeling assumptions (CORDIC without range expansion and ONE-SA's softmax reconstructed). No new physical entities are introduced.

free parameters (4)
  • tanh Chebyshev polynomial coefficients (degree 9, Q3.12, clipping [-2,2]) = not reported in paper
    Chosen by grid search over degree, bit width, and clipping range to minimize mean AE on [-4,4]; the abstract error claims depend on these coefficients.
  • sigmoid Chebyshev polynomial coefficients (degree 5, Q3.12, clipping [-3,3]) = not reported in paper
    Chosen by the same two-step search; the sigmoid error claims depend on these coefficients.
  • e^x Chebyshev polynomial coefficients (degree 3, Q3.12, clipping [-5,0]) = not reported in paper
    Chosen by ranking top candidates on softmax KL divergence for sample vectors; the softmax claims depend on these coefficients.
  • clipping thresholds x_min and x_max = tanh [-2,2], sigmoid [-3,3], e^x [-5,0]
    These bounds trade approximation error inside the interval against clipping error outside it and are chosen by search, not derived.
assumptions (5)
  • standard math Chebyshev interpolation of a continuous function gives a uniformly low-error polynomial approximation.
    Used in Sec. II to justify approximating tanh, sigmoid, and e^x by polynomials.
  • standard math Shifting softmax inputs by the maximum value leaves the output distribution unchanged.
    Used in Sec. III-C Eq. (5) to restrict e^x inputs to [-5,0].
  • domain assumption Software quantization-aware simulation accurately predicts fixed-point hardware behavior.
    The paper rounds after every operation in floating-point simulation (Sec. IV-A) and assumes this matches the synthesized 16-bit hardware.
  • ad hoc to paper The assumed ONE-SA softmax implementation is representative.
    Table II's ONE-SA comparison assumes a 4-step softmax with piecewise-linear 1/x because ONE-SA does not describe its softmax.
  • ad hoc to paper The CORDIC baseline without range expansion is a fair comparison target.
    The baseline uses CORDIC's native convergence range; the paper cites range expansion in [16] but does not apply it, which inflates the reported error gap.

how reviews work

0 comments
Cite this review

Pith. "Pith review of A Systolic Array Architecture for Nonlinear Activation Functions and Softmax Computation using Chebyshev Polynomials." pith.science (2026). https://pith.science/paper/OFU7C6C2

@misc{pith2026260804734,
  author       = {Pith},
  title        = {Pith review of: A Systolic Array Architecture for Nonlinear Activation Functions and Softmax Computation using Chebyshev Polynomials},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/OFU7C6C2}},
  note         = {Machine review of arXiv:2608.04734}
}
read the original abstract

Neural Network Accelerators have gained popularity in recent years due to their greater efficiency than CPU-based platforms. Often, these accelerators utilize different hardware units for univariate activation functions, such as tanh, and the multivariate softmax, thereby missing opportunities for resource sharing between them. In this paper, we describe a novel systolic array-based activation unit architecture that supports multiple univariate activation functions as well as the softmax function. By utilizing Chebyshev polynomial approximations, our activation function unit achieves up to 71% lower mean absolute error for tanh compared to a CORDIC baseline, while using 4.6% less area and 5.1% less power. Our softmax approximation enables a 44.6% and 79.0% lower KL divergence compared to CORDIC and a piecewise-linear approximation, respectively.

Figures

Figures reproduced from arXiv: 2608.04734 by the authors.

Figure 2
Figure 2. shows the dataflow through the systolic array on a 2× 2 PE section. Each PE receives the intermediate result IRi+1 and the value x from above and multiplies them. Afterward, it adds the next polynomial coefficient ai and propagates its result IRi and x to the PE below, also forwarding ai to its right neighbor via a register. The topmost PE receives 0 as its IR input, so its output is IRn = an, which is forwarded to … view at source ↗
Figure 1
Figure 1. Comparison of tanh(x) and e x and their polynomial approximations. Horner’s method [11] is the foundation for the efficient computation of these polynomials on a systolic array. As shown in (3), it relies on restating a polynomial as repeated multiplication of the previous intermediate result by x and subsequent addition of the next polynomial coefficient ai . a3x 3 + a2x 2 + a1x + a0 = a0 + x(a1 + x(a2 + x(a3))) (3… view at source ↗
Figure 3
Figure 3. Mean AE for tanh across [−4, 4] for the tested bitwidths with Q3.12 quantization, clipping outside of [−2, 2] for (a) and outside of [−3, 3] for (b). Values scaled by 10−1 . degree bitwidth 4 8 12 16 1 2.77 0.44 0.35 0.35 2 2.77 0.44 0.36 0.36 3 2.77 0.45 0.29 0.28 4 2.77 0.45 0.29 0.28 5 2.77 0.45 0.31 0.27 6 2.77 0.45 0.31 0.27 7 2.77 0.45 0.31 0.28 8 2.77 0.45 0.31 0.28 9 2.77 0.45 0.35 0.28 10 2.77 0.45 0.35 0.2… view at source ↗
Figures from the paper (3 more)
Figure 4
Figure 4. Figure 4: Mean AE for sigmoid across [−4, 4] for the tested bitwidths with Q3.12 quantization, clipping outside of [−2, 2] for (a) and outside of [−3, 3] for (b). Values scaled by 10−1 . integer bits used in the final quantization. In this step, for both tanh and sigmoid, the to…
Figure 5
Figure 5. Figure 5: AE of fixed-point (Q3.12) activation function imple [PITH_FULL_IMAGE:figures/full_fig_p005_5.png]
Figure 6
Figure 6. Figure 6: Area, power, and per-output latency of our designs [PITH_FULL_IMAGE:figures/full_fig_p006_6.png]

Discussion (0). Sign in to comment.

Reference graph

Works this paper leans on

16 extracted references · 13 canonical work pages

  1. [16]

    Expanding the range of con- vergence of the CORDIC algorithm,

    X. Hu, R. G. Harber, and S. C. Bass, “Expanding the range of con- vergence of the CORDIC algorithm,”IEEE Transactions on computers, vol. 40, no. 1, pp. 13–21, 2002. ©2026 IEEE. Personal use of this material is permitted. Permission from IEEE must be obtained for all other uses, in any current or future media, including reprinting/republishing this materia...

  2. [1]

    In-Datacenter Performance Analysis of a Tensor Processing Unit,

    N. P. Jouppi, C. Young, N. Patil, D. Patterson, G. Agrawal, R. Bajwa, S. Bates, S. Bhatia, N. Boden, A. Borcherset al., “In-Datacenter Performance Analysis of a Tensor Processing Unit,” inProceedings of the 44th annual international symposium on computer architecture, 2017, pp. 1–12

  3. [2]

    Flex-SFU: Activation Func- tion Acceleration with Non-Uniform Piecewise Approximation,

    R. Andri, E. Reggiani, and L. Cavigelli, “Flex-SFU: Activation Func- tion Acceleration with Non-Uniform Piecewise Approximation,”IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 2025

  4. [3]

    Design Space Exploration of Neural Network Activation Function Circuits,

    T. Yang, Y . Wei, Z. Tu, H. Zeng, M. A. Kinsy, N. Zheng, and P. Ren, “Design Space Exploration of Neural Network Activation Function Circuits,”IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, vol. 38, no. 10, pp. 1974–1978, 2018

  5. [4]

    Softermax: Hardware/Software Co-Design of an Efficient Softmax for Transformers,

    J. R. Stevens, R. Venkatesan, S. Dai, B. Khailany, and A. Raghunathan, “Softermax: Hardware/Software Co-Design of an Efficient Softmax for Transformers,” in2021 58th ACM/IEEE Design Automation Conference (DAC). IEEE, 2021, pp. 469–474

  6. [5]

    Efficient Softmax Approximation for Deep Neural Networks with Attention Mechanism,

    I. Vasyltsov and W. Chang, “Efficient Softmax Approximation for Deep Neural Networks with Attention Mechanism,”arXiv preprint arXiv:2111.10770, 2021

  7. [6]

    Aggressive Approximation of the SoftMax Function for Power-Efficient Hardware Implementations,

    F. Spagnolo, S. Perri, and P. Corsonello, “Aggressive Approximation of the SoftMax Function for Power-Efficient Hardware Implementations,” IEEE Transactions on Circuits and Systems II: Express Briefs, vol. 69, no. 3, pp. 1652–1656, 2021

  8. [7]

    Hardware-efficient Softmax Approximation for Self-Attention Networks,

    N. A. Koca, A. T. Do, and C.-H. Chang, “Hardware-efficient Softmax Approximation for Self-Attention Networks,” in2023 IEEE Interna- tional Symposium on Circuits and Systems (ISCAS). IEEE, 2023, pp. 1–5

Show all 16 references
  1. [8]

    BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding,

    J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “BERT: Pre-training of Deep Bidirectional Transformers for Language Understanding,” 2019. [Online]. Available: https://arxiv.org/abs/1810.04805

  2. [9]

    Digital Implementation of the Softmax Activation Function and the Inverse Softmax Function,

    R. S, A. S. Bharadwaj, D. S K, M. S. Khadabadi, and A. Jayaprakash, “Digital Implementation of the Softmax Activation Function and the Inverse Softmax Function,” in2022 4th International Conference on Circuits, Control, Communication and Computing (I4C), 2022, pp. 64– 67

  3. [10]

    ONE-SA: Enabling Nonlinear Operations in Systolic Arrays For Efficient and Flexible Neural Network Inference,

    R. Sun, Y . Ni, X. He, J. Zhao, and A. Zou, “ONE-SA: Enabling Nonlinear Operations in Systolic Arrays For Efficient and Flexible Neural Network Inference,”arXiv preprint arXiv:2402.00395, 2024

  4. [11]

    XXI. A new method of solving numerical equations of all orders, by continuous approximation,

    W. G. Horner, “XXI. A new method of solving numerical equations of all orders, by continuous approximation,”Philosophical Transactions of the Royal Society of London, no. 109, pp. 308–335, 1819

  5. [12]

    P. L. Chebyshev,Th ´eorie des m ´ecanismes connus sous le nom de parall´elogrammes. Imprimerie de l’Acad ´emie imp ´eriale des sciences, 1853

  6. [13]

    Goodfellow, Y

    I. Goodfellow, Y . Bengio, and A. Courville,Deep Learning. MIT Press, 2016, http://www.deeplearningbook.org

  7. [14]

    On Information and Sufficiency,

    S. Kullback and R. A. Leibler, “On Information and Sufficiency,”The annals of mathematical statistics, vol. 22, no. 1, pp. 79–86, 1951

  8. [15]

    FreePDK v2.0: Transitioning VLSI education towards nanometer variation-aware designs,

    J. E. Stine, J. Chen, I. Castellanos, G. Sundararajan, M. Qayam, P. Kumar, J. Remington, and S. Sohoni, “FreePDK v2.0: Transitioning VLSI education towards nanometer variation-aware designs,” in2009 IEEE International Conference on Microelectronic Systems Education. IEEE, 2009...

Pith tools

Reviewed August 6, 2026 · model on record in the stance chip above.