Softplus kernel

CHANGELOG.md

@@ -1,5 +1,10 @@
 # Change Log
 
+- 4.1.0
+    - Added `softplus()` and `erf()` special functions
+    - Added `sinh()` and `cosh()` trigonometric functions
+    - Added reciprocal square root `rsqrt()` function
+
 - 4.0.0
     - Data are now backed by a contiguous C buffer instead of PHP array
     - Added `fromArray()` factory method to build from PHP array

benchmarks/Unary/SoftplusMatrixBench.php

@@ -0,0 +1,62 @@
+<?php
+
+namespace Tensor\Benchmarks\Unary;
+
+use Tensor\Matrix;
+use Generator;
+
+/**
+ * Softplus as a fused kernel, against the exp/add/log composition it replaces.
+ *
+ * @Groups({"Functions"})
+ * @BeforeMethods({"setUp"})
+ */
+class SoftplusMatrixBench
+{
+    /**
+     * @var Matrix
+     */
+    protected $a;
+
+    /**
+     * @param array $params
+     */
+    public function setUp(array $params) : void
+    {
+        $this->a = Matrix::uniform(...$params['size']);
+    }
+
+    /**
+     * @return list<array{size: list<int>}>
+     */
+    public function sizes() : Generator
+    {
+        yield 'small' => ['size' => [1024, 1024]];
+        yield 'medium' => ['size' => [4096, 4096]];
+        yield 'large' => ['size' => [8192, 8192]];
+    }
+
+    /**
+     * @Subject
+     * @Iterations(5)
+     * @ParamProviders({"sizes"})
+     * @OutputTimeUnit("milliseconds", precision=3)
+     */
+    public function benchFused() : void
+    {
+        $this->a->softplus();
+    }
+
+    /**
+     * @Subject
+     * @Iterations(5)
+     * @ParamProviders({"sizes"})
+     * @OutputTimeUnit("milliseconds", precision=3)
+     */
+    public function benchComposed() : void
+    {
+        $e = $this->a->exp();
+
+        $e->add(1.0)->log();
+    }
+}

docs/Matrix.md

@@ -419,12 +419,12 @@ See [Unary](interfaces/unary.md). Each method returns a new `Matrix`.
 - `abs() : self` — absolute value of each element
 - `square() : self` — square of the matrix element-wise
 - `sqrt() : self` — square root of the matrix
+- `rsqrt() : self` — reciprocal square root of the matrix
 - `reciprocal() : self` — element-wise reciprocal of the matrix
 - `exp() : self` — exponential of the matrix
 - `expm1() : self` — exponential of each element minus 1
 - `log(float $base = M_E) : self` — logarithm of the matrix in the specified base
 - `log1p() : self` — log of 1 plus each element
-- `sigmoid() : self` — element-wise logistic function, `1 / (1 + exp(-x))`
 - `round(int $precision = 0) : self` — round the elements to a given decimal place (throws `InvalidArgumentException` if `$precision < 0`)
 - `floor() : self` — round down to the nearest integer
 - `ceil() : self` — round up to the nearest integer
@@ -445,6 +445,8 @@ See [Trigonometric](interfaces/trigonometric.md). Each method returns a new `Mat
 - `tan()` — tangent of the matrix
 - `atan()` — arc tangent of the matrix
 - `tanh() : self` — element-wise hyperbolic tangent
+- `sinh() : self` — element-wise hyperbolic sine
+- `cosh() : self` — element-wise hyperbolic cosine
 - `rad2deg()` — convert angles from radians to degrees
 - `deg2rad()` — convert angles from degrees to radians
 
@@ -466,9 +468,12 @@ See [Reductions](interfaces/reductions.md). For a `Matrix`, row-wise reductions
 
 ## Special
 
-See [Special](interfaces/special.md). These preserve the shape of the matrix without being element-wise.
+See [Special](interfaces/special.md). `softmax()` preserves the shape of the matrix; `erf()` is element-wise.
 
+- `sigmoid() : self` — element-wise logistic function, `1 / (1 + exp(-x))`
+- `softplus() : self` — element-wise softplus, `log(1 + exp(x))`
 - `softmax() : self` — exponentials normalized so each **row** sums to `1.0`; transpose on either side to normalize each column
+- `erf() : self` — element-wise Gaussian error function
 
 ## Array Access
 

docs/Vector.md

@@ -302,12 +302,14 @@ See [Unary](interfaces/unary.md). Each method returns a new `Vector`.
 - `abs() : self` — absolute value of the vector
 - `square() : self` — square the vector
 - `sqrt() : self` — square root of the vector
+- `rsqrt() : self` — reciprocal square root of the vector
 - `reciprocal() : self` — element-wise reciprocal of the vector
 - `exp() : self` — exponentiate each element
 - `expm1() : self` — exponential of each element minus 1
 - `log(float $base = M_E) : self` — log to the given base of each element
 - `log1p() : self` — log of 1 plus each element
 - `sigmoid() : self` — element-wise logistic function, `1 / (1 + exp(-x))`
+- `softplus() : self` — element-wise softplus, `log(1 + exp(x))`
 - `round(int $precision = 0) : self` — round the elements to a given decimal place (throws `InvalidArgumentException` if `$precision < 0`)
 - `floor() : self` — round down to the nearest integer
 - `ceil() : self` — round up to the nearest integer
@@ -328,6 +330,8 @@ See [Trigonometric](interfaces/trigonometric.md). Each method returns a new `Vec
 - `tan()` — tangent of the vector
 - `atan()` — arc tangent of the vector
 - `tanh() : self` — element-wise hyperbolic tangent
+- `sinh() : self` — element-wise hyperbolic sine
+- `cosh() : self` — element-wise hyperbolic cosine
 - `rad2deg()` — convert angles from radians to degrees
 - `deg2rad()` — convert angles from degrees to radians
 
@@ -348,9 +352,10 @@ See [Reductions](interfaces/reductions.md). For a `Vector` these return scalar `
 
 ## Special
 
-See [Special](interfaces/special.md). These preserve the shape of the vector without being element-wise.
+See [Special](interfaces/special.md). `softmax()` preserves the shape of the vector; `erf()` is element-wise.
 
 - `softmax() : self` — the exponentials of the vector normalized so they sum to `1.0`
+- `erf() : self` — element-wise Gaussian error function
 
 ## Array Access
 

docs/index.md

@@ -55,7 +55,7 @@ Tensor \
 | [Unary](interfaces/unary.md) | Element-wise unary functions, including clipping. |
 | [Trigonometric](interfaces/trigonometric.md) | Element-wise trigonometric functions. |
 | [Reductions](interfaces/reductions.md) | Aggregate reduction operations, including the statistical summaries. |
-| [Special](interfaces/special.md) | Operations that preserve shape without being element-wise, such as `softmax`. |
+| [Special](interfaces/special.md) | Operations that fit no other group, such as `softmax` and the error-function `erf`. |
 
 ## Decompositions
 

docs/interfaces/special.md

@@ -6,16 +6,31 @@ Operations that do not belong to any of the other interface groups.
 
 ## Overview
 
-`Special` holds the operations that neither pointwise element-wise functions (`Arithmetic`, `Comparable`, `Unary`, `Trigonometric`), nor aggregate reductions (`Reductions`), nor array-like access (`ArrayLike`) can describe.
+`Special` holds the operations that neither the basic algebraic and trigonometric families (`Arithmetic`, `Comparable`, `Unary`, `Trigonometric`), aggregate reductions (`Reductions`), nor array-like access (`ArrayLike`) can describe.
 
-`softmax()` is the sole member: it preserves the shape of the tensor, but each element of the output depends on a whole reduction axis rather than on one input element, so it is neither element-wise nor a reduction.
+`softmax()` is the only shape-preserving member whose output depends on a whole reduction axis rather than on one input element, so it is neither element-wise nor a reduction. `erf()` is an element-wise error-function transform grouped here so the error-function family stays together, distinct from the basic algebraic operations.
 
 ```php
 interface Special
 ```
 
 ## Methods
 
+### `sigmoid() : mixed`
+
+Return the element-wise logistic function, `1 / (1 + exp(-x))`.
+
+A large negative input saturates to exactly `0.0` and a large positive one to
+exactly `1.0`, so the result never overflows to an infinity.
+
+### `softplus() : mixed`
+
+Return the element-wise softplus, `log(1 + exp(x))`.
+
+The computation is numerically stable: for a large positive input it reduces to
+`x + log1p(0) == x` and for a large negative input to `log1p(0) == 0.0`, so the
+result never overflows to an infinity.
+
 ### `softmax() : mixed`
 
 Return the softmax of the tensor, i.e. the exponentials of its elements
@@ -36,3 +51,10 @@ The maximum of each row is subtracted before exponentiating, so no intermediate
 can overflow regardless of the magnitude of the input, and a single-element row
 normalizes to `1.0`. Note the degenerate case this implies: a matrix with a
 single *column* has one-element rows, so it normalizes to all ones.
+
+### `erf() : mixed`
+
+Return the element-wise Gaussian error function.
+
+A large input saturates to exactly `1.0` or `-1.0` rather than returning an
+infinity.

docs/interfaces/tensor.md

@@ -20,16 +20,15 @@ interface Tensor extends ArrayLike, Arithmetic, Comparable,
 | [ArrayLike](arraylike.md) | `shape`, `shapeString`, `size`, `map`, `reduce`, `asArray` |
 | [Arithmetic](arithmetic.md) | `multiply`, `divide`, `add`, `subtract`, `pow`, `mod` |
 | [Comparable](comparable.md) | `equal`, `notEqual`, `greater`, `greaterEqual`, `less`, `lessEqual` |
-| [Unary](unary.md) | `abs`, `square`, `sqrt`, `exp`, `expm1`, `log`, `log1p`, `sigmoid`, `round`, `floor`, `ceil`, `sign`, `negate`, `clip`, `clipLower`, `clipUpper` |
-| [Trigonometric](trigonometric.md) | `sin`, `asin`, `cos`, `acos`, `tan`, `atan`, `tanh`, `rad2deg`, `deg2rad` |
+| [Unary](unary.md) | `abs`, `square`, `sqrt`, `rsqrt`, `exp`, `expm1`, `log`, `log1p`, `round`, `floor`, `ceil`, `sign`, `negate`, `clip`, `clipLower`, `clipUpper` |
+| [Trigonometric](trigonometric.md) | `sin`, `asin`, `cos`, `acos`, `tan`, `atan`, `tanh`, `sinh`, `cosh`, `rad2deg`, `deg2rad` |
 | [Reductions](reductions.md) | `sum`, `product`, `min`, `max`, `argmin`, `argmax`, `mean`, `variance`, `median`, `quantile` |
-| [Special](special.md) | `softmax` |
+| [Special](special.md) | `sigmoid`, `softplus`, `softmax`, `erf` |
 
 Every group is homogeneous in what it does to the tensor's shape: `Arithmetic`,
 `Comparable`, `Unary`, and `Trigonometric` are element-wise and preserve it,
 `Reductions` collapse it, and `ArrayLike` is about access rather than arithmetic.
-`Special` is the escape hatch for operations that preserve the shape without being
-element-wise.
+`Special` is the escape hatch for operations that fit no other group, such as shape-preserving cross-axis operations like `softmax()` and the element-wise error-function family `erf()`.
 
 Additionally, because `ArrayLike` extends `ArrayAccess`, `IteratorAggregate`, and `Countable`, every tensor is array-accessible, iterable, and countable.
 

docs/interfaces/trigonometric.md

@@ -52,3 +52,15 @@ Return the element-wise hyperbolic tangent.
 
 A large input saturates to exactly `1.0` or `-1.0` rather than returning an
 infinity.
+
+### `sinh() : mixed`
+
+Return the element-wise hyperbolic sine.
+
+### `cosh() : mixed`
+
+Return the element-wise hyperbolic cosine.
+
+The result is always positive, and `cosh(0) == 1.0`. For a large input it grows
+as `exp(|x|) / 2`, so it can overflow to an infinity for sufficiently large
+arguments.

docs/interfaces/unary.md

@@ -26,6 +26,10 @@ Square the tensor.
 
 Return the square root of the tensor.
 
+### `rsqrt() : mixed`
+
+Return the reciprocal square root of the tensor.
+
 ### `reciprocal() : mixed`
 
 Return the element-wise reciprocal of the tensor.
@@ -49,13 +53,6 @@ Return the logarithm of the tensor in a specified base.
 
 Return the log of 1 plus the tensor, i.e. a transform.
 
-### `sigmoid() : mixed`
-
-Return the element-wise logistic function, `1 / (1 + exp(-x))`.
-
-A large negative input saturates to exactly `0.0` and a large positive one to
-exactly `1.0`, so the result never overflows to an infinity.
-
 ### `round(int $precision = 0) : self`
 
 Round the elements in the tensor to a given decimal place.

ext/include/unary.c

@@ -86,6 +86,7 @@
 
 TENSOR_UNARY_DISPATCH(abs, fabs(va[i]))
 TENSOR_UNARY_DISPATCH(sqrt, sqrt(va[i]))
+TENSOR_UNARY_DISPATCH(rsqrt, 1.0 / sqrt(va[i]))
 TENSOR_UNARY_DISPATCH(floor, floor(va[i]))
 TENSOR_UNARY_DISPATCH(ceil, ceil(va[i]))
 TENSOR_UNARY_DISPATCH(negate, -va[i])
@@ -134,13 +135,17 @@ TENSOR_UNARY(expm1, expm1(va[i]))
 TENSOR_UNARY(log, log(va[i]))
 TENSOR_UNARY(log1p, log1p(va[i]))
 TENSOR_UNARY(sigmoid, 1.0 / (1.0 + exp(-va[i])))
+TENSOR_UNARY(softplus, va[i] > 0.0 ? va[i] + log1p(exp(-va[i])) : log1p(exp(va[i])))
 TENSOR_UNARY(sin, sin(va[i]))
 TENSOR_UNARY(asin, asin(va[i]))
 TENSOR_UNARY(cos, cos(va[i]))
 TENSOR_UNARY(acos, acos(va[i]))
 TENSOR_UNARY(tan, tan(va[i]))
 TENSOR_UNARY(atan, atan(va[i]))
 TENSOR_UNARY(tanh, tanh(va[i]))
+TENSOR_UNARY(sinh, sinh(va[i]))
+TENSOR_UNARY(cosh, cosh(va[i]))
+TENSOR_UNARY(erf, erf(va[i]))
 
 #undef TENSOR_UNARY
 
@@ -462,6 +467,7 @@ void tensor_unary_dispatch_sse_init(void)
 {
 	tensor_abs_route = tensor_abs_sse;
 	tensor_sqrt_route = tensor_sqrt_sse;
+	tensor_rsqrt_route = tensor_rsqrt_sse;
 	tensor_floor_route = tensor_floor_sse;
 	tensor_ceil_route = tensor_ceil_sse;
 	tensor_negate_route = tensor_negate_sse;
@@ -480,6 +486,7 @@ void tensor_unary_dispatch_avx_init(void)
 {
 	tensor_abs_route = tensor_abs_avx;
 	tensor_sqrt_route = tensor_sqrt_avx;
+	tensor_rsqrt_route = tensor_rsqrt_avx;
 	tensor_floor_route = tensor_floor_avx;
 	tensor_ceil_route = tensor_ceil_avx;
 	tensor_negate_route = tensor_negate_avx;
@@ -498,6 +505,7 @@ void tensor_unary_dispatch_avx512_init(void)
 {
 	tensor_abs_route = tensor_abs_avx512;
 	tensor_sqrt_route = tensor_sqrt_avx512;
+	tensor_rsqrt_route = tensor_rsqrt_avx512;
 	tensor_floor_route = tensor_floor_avx512;
 	tensor_ceil_route = tensor_ceil_avx512;
 	tensor_negate_route = tensor_negate_avx512;

ext/include/unary.h

@@ -5,13 +5,18 @@
 
 void tensor_abs(zval * return_value, zval * a);
 void tensor_sqrt(zval * return_value, zval * a);
+void tensor_rsqrt(zval * return_value, zval * a);
 void tensor_exp(zval * return_value, zval * a);
 void tensor_expm1(zval * return_value, zval * a);
 void tensor_log(zval * return_value, zval * a);
 void tensor_log_base(zval * return_value, zval * a, zval * b);
 void tensor_log1p(zval * return_value, zval * a);
 void tensor_sigmoid(zval * return_value, zval * a);
+void tensor_softplus(zval * return_value, zval * a);
 void tensor_tanh(zval * return_value, zval * a);
+void tensor_sinh(zval * return_value, zval * a);
+void tensor_cosh(zval * return_value, zval * a);
+void tensor_erf(zval * return_value, zval * a);
 void tensor_sin(zval * return_value, zval * a);
 void tensor_asin(zval * return_value, zval * a);
 void tensor_cos(zval * return_value, zval * a);