Approximating Hyperbolic Tangent

The hyperbolic tangent function, ( tanh ), maps any real number to the range (-1, 1) with a smooth S-shaped curve. This property is useful as an activation function in neural networks, where it introduces non-linearity while keeping outputs bounded, and in audio signal processing, where it provides natural-sounding soft clipping for saturation and distortion effects.

In both contexts, speed matters. Neural network inference may evaluate ( tanh ) millions of times per forward pass, and audio processing demands real-time performance at sample rates of 44.1 kHz or higher. The accuracy provided by standard library implementations requires more computation than a tailored approximation.

This post surveys several approaches for approximating ( tanh ): traditional polynomial methods like Taylor series, Padé approximants, and splines, as well as more exotic techniques that exploit the IEEE-754 floating-point representation to achieve impressive speed without too much work.

§ Survey of Approximations

First, let's take a tour of some of the more common approaches to implementing fast approximations of functions like ( tanh ):

§ Taylor series

As you might recall from a distant Calculus class, a Taylor series is a polynomial expansion of a function that uses successive derivatives to generate a (infinite) summed sequence of fractions with ( x ) to successive degrees. By taking the first few terms of the sequence, we can find a rough approximation with relatively few operations.

pub fn tanhf(x: f32) -> f32 {
// Snap to 1 when polynomial begins to deviate at the tails
if x.abs() > 1.365 {
return 1f32.copysign(x);
}
let t1 = x;
let t2 = x.powi(3) * (1. / 3.);
let t3 = x.powi(5) * (2. / 15.);
let t4 = x.powi(7) * (17. / 315.);
let t5 = x.powi(9) * (62. / 2835.);
let t6 = x.powi(11) * (1382. / 155925.);
t1 - t2 + t3 - t4 + t5 - t6
}

§ Padé Approximate

Similar to the Taylor series, a Padé approximant is a ratio of two polynomials (i.e. one polynomial divided by another), which results in more accuracy, but requires more operations (including a division).

Below is a Rust adaptation of the tanh approximation in JUCE's FastMathApproximations. In mathematical parlance, this appears to be a "[7/6] Padé approximant," with a 7th-degree numerator and a 6th-degree denominator.

pub fn tanhf(x: f32) -> f32 {
// This approximation works on a limited range.
// Use input values only between -5 and +5 for limiting the error.
if x.abs() > 5. {
return 1f32.copysign(x);
}
let x2 = x * x;
let numerator = x * (135135. + x2 * (17325. + x2 * (378. + x2)));
let denominator = 135135. + x2 * (62370. + x2 * (3150. + 28. * x2));
numerator / denominator
}

§ Splines

Next, we have splines: piece-wise polynomial functions, which can be used to approximate a function by splitting it into many pieces (or subintervals) and calculating the coefficients for a (cubic) polynomial for each piece.

The example below is borrowed from Efficiently inaccurate approximation of hyperbolic tangent used as transfer function in artificial neural networks (Simos, Tsitouras), which examines the use of splines in approximating ( tanh ) for use as a transfer function in neural networks. The paper mentions that for its use, speed is the main concern, and accuracy is a lesser issue.

In tanhf3, the range [0, 18] is divided into three subintervals, and each is given its own third-degree polynomial. The specific coefficients are provided in the paper and were likely found using software like MATLAB.

pub fn tanhf3(xin: f32) -> f32 {
const N1: f32 = 0.371025186672900;
const N2: f32 = 2.572153900248530;
const N3: f32 = 18.;
match xin.abs() {
x if x <= N1 => {
-3.695076086125492e-1 * x.powi(3)
+ 1.987219343897867e-2 * x.powi(2)
+ x
}
x if x <= N2 => {
let n = x - N1;
5.928356367224758e-2 * n.powi(3)
- 3.914176949486042e-1 * n.powi(2)
+ 8.621472609449146e-1 * n
+ 3.548881072496229e-1
}
x if x <= N3 => {
let n = x - N2;
-3.347599023061577e-6 * n.powi(3)
+ 5.456777761558641e-5 * n.powi(2)
+ 7.066442941005233e-4 * n
+ 9.884026213740197e-1
}
_ => 1.,
}.copysign(xin)
}

§ Approximating By Format

Now that we've surveyed a few mathematical approaches to approximation, let's take a look at how we can use the established format of floating-point numbers (IEEE-754) to approximate ( tanh ).

Below is a diagram of the binary representation of a 32-bit floating-point number from Wikipedia. For the value of 0.15625, we can see its individual components: sign, exponent, and fraction (also known as 'mantissa').

With the binary form, the nominal value can be derived with the following equation:

$$ (s, E, M) = (−1)^s · 2^E · (1 + M/2^p), $$

where ( s, E, M, p ) are sign, (bias-added) exponent, mantissa, and the number of mantissa bits—all of which are non-negative integers.

§ K-TanH

K-TanH: Efficient TanH For Deep Learning proposed a hardware-efficient algorithm for approximating ( tanh ) using only integer operations and a small (512-bit) lookup table.

The paper provides the following pseudocode for the algorithm in Algorithm 1:

• Input:

• Input ( x_i = (s_i,E_i,M_i) )

• Parameter Tables ( T_E,T_r,T_b )

• Output:

• Output ( y_o = (s_o,E_o,M_o) )

• If ( |x_i| < T_1, y_o <- x_i )

• i.e., ( (s_o, E_o, M_o) = (s_i,E_i,M_i) )

• Else If ( |x_i| > T_2, y_o <- s_i \cdot 1 )

• i.e. ( (s_o, E_o, M_o) = (s_i,E_{bias},0) )

• Else,

• Form bit string ( t ) using lower bits of ( E_i ) and higher bits of ( M_i ).

• Fetch parameters ( \theta_t = (E_t,r_t,b_t) ) from ( T_E, T_r, T_b ) using index ( t ).

• ( s_o \leftarrow s_i, E_o \leftarrow E_t, M_o \leftarrow (M_i \gg r_t) + b_t )

• Return ( y_o )

Transforming the pseudocode above into Rust gives us the following:

pub fn tanhf(x: f32) -> f32 {
const T1: f32 = 0.25;
const T2: f32 = 3.75;
// `tanh` is symmetric around zero
let xa = x.abs();
if xa < T1 {
x
} else if xa > T2 {
1f32.copysign(x)
} else {
let x: u32 = x.to_bits();
let mi = (x >> 16) & 0b0111_1111;
let so = x & 0x8000_0000;
// t = bit string using lower bits of `E` and higher bits of `M`
let t = (x >> 20) & 0b11_111;
let (et, rt, bt) = unpack(LUT[t as usize]);
let eo = (et as u32) << 23;
let mo = (((mi >> (rt as u32)) as i32 + bt as i32) as u32) << 16;
f32::from_bits(so | eo | mo)
}
}