Skip to content

feat: add einsum() with deterministic accumulation and generic contraction kernels #34

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
CodeWithKyrian merged 1 commit into from
Jun 15, 2026
Merged

feat: add einsum() with deterministic accumulation and generic contraction kernels #34

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
43 changes: 43 additions & 0 deletions SPEC.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -362,6 +362,49 @@ Sort kind selection is enum-based via `SortKind`:
- [x] 9.4.1: `$a->solve($b)` - Solve linear equations
- [x] 9.4.2: `$a->lstsq($b)` - Least squares solution

### 9.5 Einstein Summation (REQ-9.5)
**Priority**: LOW

Einstein summation notation for two-operand tensor operations with deterministic accumulation order (no BLAS tiling).

**Subscript syntax:** `"labels_left->labels_right"` where each side is comma-separated lowercase letters. When `->` is omitted, output labels are inferred from labels appearing exactly once across inputs.

**Shape semantics:** Shared labels across inputs must match dimensions. Labels appearing only on output are broadcast axes (size 1). Labels appearing twice without appearing in output are summed over (contracted).

**`optimize` parameter:** Not applicable. All contractions use fixed canonical loop order (output axes outermost left-to-right, contracted axes innermost).

**Supported patterns:**

| Pattern | Subscripts | Description |
|---|---|---|
| Dot product | `i,i->` | Two 1D vectors → scalar |
| Element-wise | `i,i->i`, `ij,ij->ij` | Same-shape element-wise |
| Outer product | `i,j->ij` | 1D × 1D → 2D |
| Matrix-vector | `ij,j->i` | 2D × 1D → 1D |
| Vector-matrix | `i,ij->j` | 1D × 2D → 1D |
| Matrix multiply | `ij,jk->ik` | 2D × 2D → 2D (deterministic triple-loop) |
| Matrix × transposed | `ij,kj->ik` | 2D × 2D → 2D |
| Hadamard product | `ij,ij->ij` | Element-wise multiply |
| Trace | `ii->` | Sum diagonal → scalar |
| Diagonal | `ii->i` | Extract diagonal → 1D |
| Transpose | `ij->ji` | Single-operand axis swap |
| Sum over axis | `ij->i`, `ij->j`, `i->` | Reduction |

A generic contraction engine handles any pattern not covered by optimized kernels.

**Requirements:**
- [x] 9.5.1: `$a->einsum($subscripts, $b)` — Two-operand einsum
- [x] 9.5.2: Implicit output mode when `->` is omitted
- [x] 9.5.3: Shape validation — shared labels must match dimensions
- [x] 9.5.4: Label validation — no label may appear 3+ times across inputs
- [x] 9.5.5: DType promotion using existing `promote()` rules
- [x] 9.5.6: Deterministic accumulation — fixed loop order, no BLAS
- [x] 9.5.7: Optimized kernel for `ij,jk->ik` (gemm_ordered)
- [x] 9.5.8: Optimized kernel for `i,i->` (dot_ordered)
- [x] 9.5.9: Optimized kernel for `i,j->ij` (outer_ordered)
- [x] 9.5.10: Generic contraction engine for unoptimized patterns
- [x] 9.5.11: Global function alias `einsum()` in Functions.php

## 10. Iteration and Application

### 10.1 Iterator Support (REQ-10.1)
Expand Down
5 changes: 4 additions & 1 deletion docs/api/global-functions.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,6 +48,7 @@ You can combine imports from the base namespace and sub-namespaces in one `use f
| `full` | `NDArray::full()` | [Array Creation — full](/api/array-creation#ndarray-full) |
| `from_buffer` | `NDArray::fromBuffer()` | [Array Creation — fromBuffer](/api/array-creation#ndarray-frombuffer) |
| `from_bytes` | `NDArray::fromBytes()` | [Array Creation — fromBytes](/api/array-creation#ndarray-frombytes) |
| `from_scalar` | `NDArray::fromScalar()` | [Array Creation — fromScalar](/api/array-creation#ndarray-fromscalar) |
| `zeros_like` | `NDArray::zerosLike()` | [Array Creation](/api/array-creation) |
| `ones_like` | `NDArray::onesLike()` | [Array Creation](/api/array-creation) |
| `full_like` | `NDArray::fullLike()` | [Array Creation](/api/array-creation) |
Expand All @@ -65,7 +66,8 @@ You can combine imports from the base namespace and sub-namespaces in one `use f
| `tile` | `NDArray::tile()` | [Array Manipulation - tile](/api/array-manipulation#tile) |
| `repeat` | `NDArray::repeat()` | [Array Manipulation - repeat](/api/array-manipulation#repeat) |
| `copy` | `$a->copy()` | [Array Manipulation](/api/array-manipulation) |
| `astype` | `$a->astype()` | [Array Manipulation](/api/array-manipulation) |
| `astype` | `$a->astype()` | [Array Manipulation](/api/array-manipulation) |
| `cast` | `$a->cast()` | [Array Manipulation](/api/array-manipulation) |

### Element-wise math and arithmetic

Expand Down Expand Up @@ -246,6 +248,7 @@ use function PhpMlKit\NDArray\Linalg\matmul;
| `norm` | `$a->norm()` | [Linear Algebra – norm](/api/linear-algebra#norm) |
| `dot` | `$a->dot()` | [Linear Algebra – dot](/api/linear-algebra#dot) |
| `matmul` | `$a->matmul()` | [Linear Algebra – matmul](/api/linear-algebra#matmul) |
| `einsum` | `$a->einsum()` | [Linear Algebra – einsum](/api/linear-algebra#einsum) |
| `diagonal` | `$a->diagonal()` | [Linear Algebra – diagonal](/api/linear-algebra#diagonal) |
| `diag` | `$a->diag()` | [Linear Algebra – diag](/api/linear-algebra#diag) |
| `trace` | `$a->trace()` | [Linear Algebra – trace](/api/linear-algebra#trace) |
Expand Down
66 changes: 66 additions & 0 deletions docs/api/linear-algebra.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -167,6 +167,72 @@ print_r($c->toArray());
// Output: [[19, 22], [43, 50]]
```

## einsum()

```php
public function einsum(string $subscripts, ?NDArray $other = null): NDArray
```

Einstein summation with deterministic accumulation order. Evaluates the subscript expression using fixed nested loops — no BLAS tiling, identical results on every call.

::: info Deterministic accumulation
einsum does not use BLAS for any pattern, including matrix multiplication. Every loop iterates in a fixed canonical order, guaranteeing bit-identical output on every run. For the same operation with BLAS-speed (but non-deterministic accumulation), use the equivalent dedicated method instead — e.g. `matmul()` for `ij,jk->ik`, `dot()` for `i,i->`, or `multiply()` for element-wise `ij,ij->ij`.
:::

**Supported subscript patterns:**

| Pattern | Subscripts | Description |
|---|---|---|
| Matrix multiply | `ij,jk->ik` | 2D × 2D → 2D |
| Matrix × transposed | `ij,kj->ik` | 2D × 2D → 2D |
| Matrix-vector | `ij,j->i` | 2D × 1D → 1D |
| Vector-matrix | `i,ij->j` | 1D × 2D → 1D |
| Dot product | `i,i->` | Two 1D vectors → scalar |
| Outer product | `i,j->ij` | 1D × 1D → 2D |
| Element-wise | `ij,ij->ij`, `i,i->i` | Same-shape element-wise |
| Trace | `ii->` | Sum diagonal → scalar |
| Diagonal | `ii->i` | Extract diagonal → 1D |
| Transpose | `ij->ji` | Single-operand axis swap |
| Sum over axis | `ij->i`, `ij->j` | Reduction over one axis |
| Sum all | `i->` | Sum all elements → scalar |

When `->` is omitted, output labels are inferred from labels appearing exactly once across all operands (e.g. `"ij,jk"` → `"ik"`). Single-operand patterns like `"ii->"` omit the second argument.

### Parameters

| Parameter | Type | Description |
|-----------|------|-------------|
| `$subscripts` | `string` | Einstein summation subscript (e.g. `"ij,jk->ik"`). |
| `$other` | `NDArray` | The second operand. |

### Returns

- `NDArray` - Result of the contraction.

### Examples

```php
$a = NDArray::array([[1, 2, 3], [4, 5, 6]]); // 2×3
$b = NDArray::array([[1, 2], [3, 4], [5, 6]]); // 3×2

// Matrix multiplication
$c = $a->einsum('ij,jk->ik', $b);
print_r($c->toArray());
// Output: [[22, 28], [49, 64]]

// Dot product
$x = NDArray::array([1, 2, 3]);
$y = NDArray::array([4, 5, 6]);
$d = $x->einsum('i,i->', $y);
echo $d->toScalar();
// Output: 32

// Outer product
$o = NDArray::array([1, 2])->einsum('i,j->ij', NDArray::array([3, 4, 5]));
print_r($o->toArray());
// Output: [[3, 4, 5], [6, 8, 10]]
```

## diagonal()

```php
Expand Down
14 changes: 14 additions & 0 deletions include/ndarray_php.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -567,6 +567,20 @@ int32_t ndarray_ne_scalar(const struct NdArrayHandle *a,
uintptr_t *out_shape,
uintptr_t max_ndim);

/**
* FFI entry point for einsum (1 or 2 operands, null b for single-op).
*/
int32_t ndarray_einsum(const struct NdArrayHandle *a,
const struct ArrayMetadata *a_meta,
const struct NdArrayHandle *b,
const struct ArrayMetadata *b_meta,
const char *subscripts,
struct NdArrayHandle **out_handle,
uint8_t *out_dtype,
uintptr_t *out_ndim,
uintptr_t *out_shape_ptr,
uintptr_t max_ndim);

/**
* One-dimensional complex FFT along `axis`. Real inputs are promoted to complex.
* `n == 0` keeps the current length along `axis` (optional padding length).
Expand Down
67 changes: 51 additions & 16 deletions rust/src/ffi/arithmetic/maximum.rs
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,8 @@ use crate::helpers::elementwise_minmax::ElementwiseMaximum;
use crate::helpers::error::{set_last_error, ERR_GENERIC, SUCCESS};
use crate::helpers::write_output_metadata;
use crate::helpers::{
extract_array_f32, extract_array_f64, extract_array_i16, extract_array_i32,
extract_array_i64, extract_array_i8, extract_array_u16, extract_array_u32,
extract_array_u64, extract_array_u8,
extract_array_f32, extract_array_f64, extract_array_i16, extract_array_i32, extract_array_i64,
extract_array_i8, extract_array_u16, extract_array_u32, extract_array_u64, extract_array_u8,
};
use crate::types::dtype::DType;
use crate::types::{ArrayData, ArrayMetadata, NDArrayWrapper, NdArrayHandle};
Expand Down Expand Up @@ -83,8 +82,12 @@ pub unsafe extern "C" fn ndarray_maximum_scalar(
out_shape: *mut usize,
max_ndim: usize,
) -> i32 {
if a.is_null() || meta.is_null() || out.is_null() || out_dtype.is_null()
|| out_ndim.is_null() || out_shape.is_null()
if a.is_null()
|| meta.is_null()
|| out.is_null()
|| out_dtype.is_null()
|| out_ndim.is_null()
|| out_shape.is_null()
{
return ERR_GENERIC;
}
Expand All @@ -100,15 +103,21 @@ pub unsafe extern "C" fn ndarray_maximum_scalar(
return ERR_GENERIC;
};
let result = arr.mapv(|x| x.max(scalar));
NDArrayWrapper { data: ArrayData::Float64(Arc::new(RwLock::new(result))), dtype: DType::Float64 }
NDArrayWrapper {
data: ArrayData::Float64(Arc::new(RwLock::new(result))),
dtype: DType::Float64,
}
}
DType::Float32 => {
let Some(arr) = extract_array_f32(a_wrapper, meta) else {
set_last_error("Failed to extract f32 array".to_string());
return ERR_GENERIC;
};
let result = arr.mapv(|x| x.max(scalar as f32));
NDArrayWrapper { data: ArrayData::Float32(Arc::new(RwLock::new(result))), dtype: DType::Float32 }
NDArrayWrapper {
data: ArrayData::Float32(Arc::new(RwLock::new(result))),
dtype: DType::Float32,
}
}
DType::Int64 => {
let Some(arr) = extract_array_i64(a_wrapper, meta) else {
Expand All @@ -117,7 +126,10 @@ pub unsafe extern "C" fn ndarray_maximum_scalar(
};
let s = scalar as i64;
let result = arr.mapv(|x| x.max(s));
NDArrayWrapper { data: ArrayData::Int64(Arc::new(RwLock::new(result))), dtype: DType::Int64 }
NDArrayWrapper {
data: ArrayData::Int64(Arc::new(RwLock::new(result))),
dtype: DType::Int64,
}
}
DType::Int32 => {
let Some(arr) = extract_array_i32(a_wrapper, meta) else {
Expand All @@ -126,7 +138,10 @@ pub unsafe extern "C" fn ndarray_maximum_scalar(
};
let s = scalar as i32;
let result = arr.mapv(|x| x.max(s));
NDArrayWrapper { data: ArrayData::Int32(Arc::new(RwLock::new(result))), dtype: DType::Int32 }
NDArrayWrapper {
data: ArrayData::Int32(Arc::new(RwLock::new(result))),
dtype: DType::Int32,
}
}
DType::Int16 => {
let Some(arr) = extract_array_i16(a_wrapper, meta) else {
Expand All @@ -135,7 +150,10 @@ pub unsafe extern "C" fn ndarray_maximum_scalar(
};
let s = scalar as i16;
let result = arr.mapv(|x| x.max(s));
NDArrayWrapper { data: ArrayData::Int16(Arc::new(RwLock::new(result))), dtype: DType::Int16 }
NDArrayWrapper {
data: ArrayData::Int16(Arc::new(RwLock::new(result))),
dtype: DType::Int16,
}
}
DType::Int8 => {
let Some(arr) = extract_array_i8(a_wrapper, meta) else {
Expand All @@ -144,7 +162,10 @@ pub unsafe extern "C" fn ndarray_maximum_scalar(
};
let s = scalar as i8;
let result = arr.mapv(|x| x.max(s));
NDArrayWrapper { data: ArrayData::Int8(Arc::new(RwLock::new(result))), dtype: DType::Int8 }
NDArrayWrapper {
data: ArrayData::Int8(Arc::new(RwLock::new(result))),
dtype: DType::Int8,
}
}
DType::Uint64 => {
let Some(arr) = extract_array_u64(a_wrapper, meta) else {
Expand All @@ -153,7 +174,10 @@ pub unsafe extern "C" fn ndarray_maximum_scalar(
};
let s = (scalar.max(0.0)) as u64;
let result = arr.mapv(|x| x.max(s));
NDArrayWrapper { data: ArrayData::Uint64(Arc::new(RwLock::new(result))), dtype: DType::Uint64 }
NDArrayWrapper {
data: ArrayData::Uint64(Arc::new(RwLock::new(result))),
dtype: DType::Uint64,
}
}
DType::Uint32 => {
let Some(arr) = extract_array_u32(a_wrapper, meta) else {
Expand All @@ -162,7 +186,10 @@ pub unsafe extern "C" fn ndarray_maximum_scalar(
};
let s = (scalar.max(0.0)) as u32;
let result = arr.mapv(|x| x.max(s));
NDArrayWrapper { data: ArrayData::Uint32(Arc::new(RwLock::new(result))), dtype: DType::Uint32 }
NDArrayWrapper {
data: ArrayData::Uint32(Arc::new(RwLock::new(result))),
dtype: DType::Uint32,
}
}
DType::Uint16 => {
let Some(arr) = extract_array_u16(a_wrapper, meta) else {
Expand All @@ -171,7 +198,10 @@ pub unsafe extern "C" fn ndarray_maximum_scalar(
};
let s = (scalar.max(0.0)) as u16;
let result = arr.mapv(|x| x.max(s));
NDArrayWrapper { data: ArrayData::Uint16(Arc::new(RwLock::new(result))), dtype: DType::Uint16 }
NDArrayWrapper {
data: ArrayData::Uint16(Arc::new(RwLock::new(result))),
dtype: DType::Uint16,
}
}
DType::Uint8 => {
let Some(arr) = extract_array_u8(a_wrapper, meta) else {
Expand All @@ -180,7 +210,10 @@ pub unsafe extern "C" fn ndarray_maximum_scalar(
};
let s = (scalar.max(0.0)) as u8;
let result = arr.mapv(|x| x.max(s));
NDArrayWrapper { data: ArrayData::Uint8(Arc::new(RwLock::new(result))), dtype: DType::Uint8 }
NDArrayWrapper {
data: ArrayData::Uint8(Arc::new(RwLock::new(result))),
dtype: DType::Uint8,
}
}
DType::Bool => {
set_last_error("maximum_scalar() not supported for Bool type".to_string());
Expand All @@ -192,7 +225,9 @@ pub unsafe extern "C" fn ndarray_maximum_scalar(
}
};

if let Err(e) = write_output_metadata(&result_wrapper, out_dtype, out_ndim, out_shape, max_ndim) {
if let Err(e) =
write_output_metadata(&result_wrapper, out_dtype, out_ndim, out_shape, max_ndim)
{
set_last_error(e);
return ERR_GENERIC;
}
Expand Down
Loading