phpmlkit · CodeWithKyrian · Apr 19, 2026 · Apr 19, 2026
diff --git a/SPEC.md b/SPEC.md
@@ -37,7 +37,7 @@ The primary array type representing N-dimensional arrays.
 - [x] 2.2.2: `uint8`, `uint16`, `uint32`, `uint64` - Unsigned integers
 - [x] 2.2.3: `float32`, `float64` - Floating-point numbers
 - [x] 2.2.4: `bool` - Boolean values
-- [ ] 2.2.5: `complex64`, `complex128` - Complex numbers (future)
+- [x] 2.2.5: `complex64`, `complex128` - Complex numbers
 
 ### 2.3 Array Properties (REQ-2.3)
 **Priority**: CRITICAL
@@ -337,8 +337,11 @@ Sort kind selection is enum-based via `SortKind`:
 - [x] 9.2.1: `$array->svd($full_matrices = true)` - Singular value decomposition
 - [x] 9.2.2: `$array->qr()` - QR decomposition
 - [x] 9.2.3: `$array->cholesky()` - Cholesky decomposition
-- [ ] 9.2.4: `$array->eig()` - Eigenvalue decomposition (blocked on complex number support)
-- [ ] 9.2.5: `$array->lu()` - LU decomposition (not exposed; used internally by solve/inv/det)
+- [x] 9.2.4: `$array->eig()` - General eigenvalue decomposition (eigenvalues and eigenvectors)
+- [x] 9.2.5: `$array->eigvals()` - Eigenvalues only (general matrix)
+- [x] 9.2.6: `$array->eigh($upper = false)` - Hermitian/symmetric eigen decomposition
+- [x] 9.2.7: `$array->eigvalsh($upper = false)` - Hermitian/symmetric eigenvalues only
+- [ ] 9.2.8: `$array->lu()` - LU decomposition (not exposed; used internally by solve/inv/det)
 
 ### 9.3 Matrix Properties (REQ-9.3)
 **Priority**: MEDIUM (requires BLAS)
@@ -423,7 +426,7 @@ Sort kind selection is enum-based via `SortKind`:
 - [x] 12.2.2: Dot products (DOT)
 - [x] 12.2.3: Matrix-vector operations (GEMV)
 - [x] 12.2.4: Linear system solving
-- [x] 12.2.5: SVD decompositions (eigenvalue blocked on complex number support)
+- [x] 12.2.5: SVD decompositions (complex support available, eigenvalue decomposition pending)
 
 ## 13. Memory Management
 
@@ -507,7 +510,7 @@ Sort kind selection is enum-based via `SortKind`:
 
 **Requirements**:
 - [ ] 16.1.1: Unit tests for all public methods
-- [ ] 16.1.2: Integration tests for complex workflows
+- [x] 16.1.2: Integration tests for complex workflows
 - [ ] 16.1.3: Property-based testing for mathematical correctness
 - [ ] 16.1.4: Benchmark suite for performance regression
 - [ ] 16.1.5: Memory leak detection tests
@@ -570,7 +573,7 @@ Sort kind selection is enum-based via `SortKind`:
 
 - [ ] 19.1.1: GPU acceleration via CUDA/ROCm
 - [ ] 19.1.2: Sparse matrix support
-- [ ] 19.1.3: Complex number support
+- [x] 19.1.3: Complex number support
 - [ ] 19.1.4: FFT operations
 - [ ] 19.1.5: Image processing utilities
 - [ ] 19.1.6: Pandas-like DataFrame functionality

diff --git a/docs/api/linear-algebra.md b/docs/api/linear-algebra.md
@@ -303,6 +303,114 @@ $a = NDArray::array([
 $reconstructed = $q->matmul($r);
 ```
 
+## eig()
+
+```php
+public function eig(): array
+```
+
+Eigenvalue decomposition of a square matrix: `A v = λ v`.
+
+For real `Float32` / `Float64` input, eigenvalues and eigenvectors are complex. For complex input, output dtypes match the array dtype.
+
+### Returns
+
+- `array{0: NDArray, 1: NDArray}` - `[eigenvalues, eigenvectors]`
+
+### Examples
+
+```php
+$a = NDArray::array([
+    [4.0, 2.0],
+    [1.0, 3.0],
+]);
+
+[$w, $v] = $a->eig();
+// $w: eigenvalues (complex dtype for real input)
+// $v: eigenvector columns
+```
+
+## eigvals()
+
+```php
+public function eigvals(): NDArray
+```
+
+Eigenvalues only for a general square matrix (no eigenvectors).
+
+For real input, eigenvalues are complex. For complex input, output dtype matches the array dtype.
+
+### Returns
+
+- `NDArray` - 1D vector of eigenvalues
+
+### Examples
+
+```php
+$a = NDArray::array([[4.0, 2.0], [1.0, 3.0]]);
+$w = $a->eigvals();
+```
+
+## eigh()
+
+```php
+public function eigh(bool $upper = false): array
+```
+
+Eigen decomposition for Hermitian (or real symmetric) matrices.
+
+Eigenvalues are real. Eigenvectors use the same element type as `A`. Only the stored triangle is read: lower if `$upper` is false, upper if true.
+
+### Parameters
+
+| Parameter | Type | Description |
+|-----------|------|-------------|
+| `$upper` | `bool` | If true, read the upper triangle; if false, the lower. Optional. Default: `false`. |
+
+### Returns
+
+- `array{0: NDArray, 1: NDArray}` - `[eigenvalues, eigenvectors]`
+
+### Examples
+
+```php
+$a = NDArray::array([
+    [2.0, 1.0],
+    [1.0, 2.0],
+]);
+
+[$w, $v] = $a->eigh();
+```
+
+## eigvalsh()
+
+```php
+public function eigvalsh(bool $upper = false): NDArray
+```
+
+Eigenvalues only for a Hermitian (or real symmetric) matrix. Same triangle convention as `eigh()`.
+
+### Parameters
+
+| Parameter | Type | Description |
+|-----------|------|-------------|
+| `$upper` | `bool` | If true, read the upper triangle; if false, the lower. Optional. Default: `false`. |
+
+### Returns
+
+- `NDArray` - 1D vector of real eigenvalues
+
+### Examples
+
+```php
+$a = NDArray::array([
+    [2.0, 1.0],
+    [1.0, 2.0],
+]);
+
+$w = $a->eigvalsh();
+```
+
 ## cholesky()
 
 ```php

diff --git a/docs/guide/fundamentals/operations.md b/docs/guide/fundamentals/operations.md
@@ -203,6 +203,10 @@ $v1->dot($v2);      // 32 (scalar)
 
 // Singular Value Decomposition
 [$u, $s, $vt] = $a->svd();
+
+// Eigen decomposition (square matrix)
+$a = NDArray::array([[4.0, 2.0], [1.0, 3.0]]);
+[$eigvals, $eigvecs] = $a->eig();
 ```
 
 **Important Distinction:**
@@ -214,7 +218,7 @@ $v1->dot($v2);      // 32 (scalar)
 - `transpose()` - Transpose matrix
 - `diagonal()`, `trace()` - Diagonal operations
 - `norm()` - Vector/matrix norms
-- `svd()`, `qr()`, `cholesky()` - Matrix decompositions
+- `svd()`, `qr()`, `eig()`, `eigvals()`, `eigh()`, `eigvalsh()`, `cholesky()` - Matrix decompositions
 - `inv()`, `det()`, `pinv()`, `cond()`, `rank()` - Matrix properties
 - `solve()`, `lstsq()` - Linear system solvers
 

diff --git a/docs/guide/getting-started/numpy-migration.md b/docs/guide/getting-started/numpy-migration.md
@@ -167,6 +167,10 @@ $mask = $a->gt($b);
 | `a.T` | `$a->transpose()` | Use method instead |
 | `np.linalg.svd(a)` | `$a->svd()` | Singular value decomposition |
 | `np.linalg.qr(a)` | `$a->qr()` | QR decomposition |
+| `np.linalg.eig(a)` | `$a->eig()` | Eigenvalues and eigenvectors |
+| `np.linalg.eigvals(a)` | `$a->eigvals()` | Eigenvalues only |
+| `np.linalg.eigh(a)` | `$a->eigh()` | Hermitian / symmetric |
+| `np.linalg.eigvalsh(a)` | `$a->eigvalsh()` | Hermitian eigenvalues only |
 | `np.linalg.cholesky(a)` | `$a->cholesky()` | Cholesky decomposition |
 | `np.linalg.inv(a)` | `$a->inv()` | Matrix inverse |
 | `np.linalg.det(a)` | `$a->det()` | Determinant |

diff --git a/include/ndarray_php.h b/include/ndarray_php.h
@@ -904,6 +904,94 @@ int32_t ndarray_dot(const struct NdArrayHandle *a,
                     uintptr_t *out_shape,
                     uintptr_t max_ndim);
 
+/**
+ * Compute eigenvalue decomposition: `A * v = lambda * v`.
+ *
+ * For real input, eigenvalues and eigenvectors are complex.
+ * For complex input, output type matches input type.
+ */
+int32_t ndarray_eig(const struct NdArrayHandle *a,
+                    const struct ArrayMetadata *a_meta,
+                    struct NdArrayHandle **out_eigvals,
+                    uint8_t *out_dtype_eigvals,
+                    uintptr_t *out_ndim_eigvals,
+                    uintptr_t *out_shape_eigvals,
+                    uintptr_t max_ndim,
+                    struct NdArrayHandle **out_eigvecs,
+                    uint8_t *out_dtype_eigvecs,
+                    uintptr_t *out_ndim_eigvecs,
+                    uintptr_t *out_shape_eigvecs);
+
+/**
+ * Compute eigenvalue decomposition for Hermitian/symmetric matrix.
+ *
+ * Eigenvalues are always real. Eigenvectors have the same type as input.
+ *
+ * # Arguments
+ * * `a` - Input matrix handle (n x n, 2D, square, Hermitian/symmetric)
+ * * `a_meta` - Array metadata
+ * * `uplo` - 0 for Lower triangle, 1 for Upper triangle
+ * * `out_eigvals` - Output eigenvalues vector handle (length n, real)
+ * * `out_dtype_eigvals` - Eigenvalues dtype output
+ * * `out_ndim_eigvals` - Eigenvalues ndim output
+ * * `out_shape_eigvals` - Eigenvalues shape output
+ * * `max_ndim` - Maximum number of dimensions
+ * * `out_eigvecs` - Output eigenvectors matrix handle (n x n)
+ * * `out_dtype_eigvecs` - Eigenvectors dtype output
+ * * `out_ndim_eigvecs` - Eigenvectors ndim output
+ * * `out_shape_eigvecs` - Eigenvectors shape output
+ */
+int32_t ndarray_eigh(const struct NdArrayHandle *a,
+                     const struct ArrayMetadata *a_meta,
+                     uint8_t uplo,
+                     struct NdArrayHandle **out_eigvals,
+                     uint8_t *out_dtype_eigvals,
+                     uintptr_t *out_ndim_eigvals,
+                     uintptr_t *out_shape_eigvals,
+                     uintptr_t max_ndim,
+                     struct NdArrayHandle **out_eigvecs,
+                     uint8_t *out_dtype_eigvecs,
+                     uintptr_t *out_ndim_eigvecs,
+                     uintptr_t *out_shape_eigvecs);
+
+/**
+ * Compute eigenvalues only (no eigenvectors) for a general matrix.
+ *
+ * For real input, eigenvalues are complex.
+ * For complex input, output type matches input type.
+ */
+int32_t ndarray_eigvals(const struct NdArrayHandle *a,
+                        const struct ArrayMetadata *a_meta,
+                        struct NdArrayHandle **out_eigvals,
+                        uint8_t *out_dtype_eigvals,
+                        uintptr_t *out_ndim_eigvals,
+                        uintptr_t *out_shape_eigvals,
+                        uintptr_t max_ndim);
+
+/**
+ * Compute eigenvalues only (no eigenvectors) for a Hermitian/symmetric matrix.
+ *
+ * Eigenvalues are always real.
+ *
+ * # Arguments
+ * * `a` - Input matrix handle (n x n, 2D, square, Hermitian/symmetric)
+ * * `a_meta` - Array metadata
+ * * `uplo` - 0 for Lower triangle, 1 for Upper triangle
+ * * `out_eigvals` - Output eigenvalues vector handle (length n, real)
+ * * `out_dtype_eigvals` - Eigenvalues dtype output
+ * * `out_ndim_eigvals` - Eigenvalues ndim output
+ * * `out_shape_eigvals` - Eigenvalues shape output
+ * * `max_ndim` - Maximum number of dimensions
+ */
+int32_t ndarray_eigvalsh(const struct NdArrayHandle *a,
+                         const struct ArrayMetadata *a_meta,
+                         uint8_t uplo,
+                         struct NdArrayHandle **out_eigvals,
+                         uint8_t *out_dtype_eigvals,
+                         uintptr_t *out_ndim_eigvals,
+                         uintptr_t *out_shape_eigvals,
+                         uintptr_t max_ndim);
+
 /**
  * Create a 2D array with the given 1D data on a diagonal.
  *