From 355ce4dd0be9f3eb484d5310dc803fa92a9a17e8 Mon Sep 17 00:00:00 2001 From: kelikovic Date: Sun, 13 Sep 2026 03:32:08 +0200 Subject: [PATCH 1/3] feat: statically-typed matrix/vector ops and explicit axis broadcasting Add scale() and axis-explicit addVector/subtractVector/multiplyVector/divideVector to Matrix64/Matrix32, plus Vector64/Vector32 scale(), giving a fully typed multiply path and removing the square-matrix broadcast ambiguity. --- lib/src/matrix.dart | 47 +++++++++++++++++++++++++++++ lib/src/matrix32.dart | 47 +++++++++++++++++++++++++++++ lib/src/vector.dart | 8 +++++ lib/src/vector32.dart | 8 +++++ test/matrix_api_test.dart | 63 +++++++++++++++++++++++++++++++++++++++ 5 files changed, 173 insertions(+) diff --git a/lib/src/matrix.dart b/lib/src/matrix.dart index 8e7883c..c475e5d 100644 --- a/lib/src/matrix.dart +++ b/lib/src/matrix.dart @@ -653,6 +653,29 @@ class Matrix64 { Matrix64 hadamard(Object other) => _binary(other, ElementwiseOp.multiply); + /// Scalar multiply that preserves the static [Matrix64] type. + /// + /// `operator *` returns `dynamic` because it also dispatches to + /// matrix/vector products; use [scale] when the result must stay + /// statically typed for chaining. + Matrix64 scale(num factor) => _binary(factor, ElementwiseOp.multiply); + + /// Add [vector] to every line along [axis] with an explicit + /// orientation. `Axis.columns` expects `vector.length == columns`, + /// `Axis.rows` expects `vector.length == rows` — removing the + /// square-matrix ambiguity of the broadcasting operators. + Matrix64 addVector(Vector64 vector, Axis axis) => + _broadcastAxis(vector, axis, ElementwiseOp.add); + + Matrix64 subtractVector(Vector64 vector, Axis axis) => + _broadcastAxis(vector, axis, ElementwiseOp.subtract); + + Matrix64 multiplyVector(Vector64 vector, Axis axis) => + _broadcastAxis(vector, axis, ElementwiseOp.multiply); + + Matrix64 divideVector(Vector64 vector, Axis axis) => + _broadcastAxis(vector, axis, ElementwiseOp.divide); + Matrix64 matmul(Matrix64 other) { if (columns != other.rows) { throw ArgumentError( @@ -941,6 +964,30 @@ class Matrix64 { ); } + Matrix64 _broadcastAxis(Vector64 vector, Axis axis, ElementwiseOp op) { + final expected = axis == Axis.columns ? columns : rows; + if (vector.length != expected) { + throw ArgumentError( + 'Vector64 length ${vector.length} must equal ' + '${axis == Axis.columns ? 'columns $columns' : 'rows $rows'} ' + 'to broadcast along $axis.', + ); + } + final vectorData = vector.unsafeValuesView; + final result = Float64List(count); + for (var row = 0; row < rows; row++) { + final offset = row * columns; + for (var column = 0; column < columns; column++) { + result[offset + column] = _apply( + _data[offset + column], + vectorData[axis == Axis.columns ? column : row], + op, + ); + } + } + return Matrix64.storage(rows, columns, result, copy: false); + } + Matrix64 _broadcastVector(Vector64 vector, ElementwiseOp op) { final vectorData = vector.values; final result = Float64List(count); diff --git a/lib/src/matrix32.dart b/lib/src/matrix32.dart index 9184db5..bb484ed 100644 --- a/lib/src/matrix32.dart +++ b/lib/src/matrix32.dart @@ -396,6 +396,29 @@ class Matrix32 { Matrix32 hadamard(Object other) => _binary(other, ElementwiseOp.multiply); + /// Scalar multiply that preserves the static [Matrix32] type. + /// + /// `operator *` returns `dynamic` because it also dispatches to + /// matrix/vector products; use [scale] when the result must stay + /// statically typed for chaining. + Matrix32 scale(num factor) => _binary(factor, ElementwiseOp.multiply); + + /// Add [vector] to every line along [axis] with an explicit + /// orientation. `Axis.columns` expects `vector.length == columns`, + /// `Axis.rows` expects `vector.length == rows` — removing the + /// square-matrix ambiguity of the broadcasting operators. + Matrix32 addVector(Vector32 vector, Axis axis) => + _broadcastAxis32(vector, axis, ElementwiseOp.add); + + Matrix32 subtractVector(Vector32 vector, Axis axis) => + _broadcastAxis32(vector, axis, ElementwiseOp.subtract); + + Matrix32 multiplyVector(Vector32 vector, Axis axis) => + _broadcastAxis32(vector, axis, ElementwiseOp.multiply); + + Matrix32 divideVector(Vector32 vector, Axis axis) => + _broadcastAxis32(vector, axis, ElementwiseOp.divide); + Matrix32 vstack(Matrix32 other) => Matrix32.fromMatrix(toFloat64().vstack(other.toFloat64())); @@ -719,6 +742,30 @@ class Matrix32 { ); } + Matrix32 _broadcastAxis32(Vector32 vector, Axis axis, ElementwiseOp op) { + final expected = axis == Axis.columns ? columns : rows; + if (vector.length != expected) { + throw ArgumentError( + 'Vector32 length ${vector.length} must equal ' + '${axis == Axis.columns ? 'columns $columns' : 'rows $rows'} ' + 'to broadcast along $axis.', + ); + } + final vectorData = vector.unsafeValuesView; + final result = Float32List(count); + for (var row = 0; row < rows; row++) { + final offset = row * columns; + for (var column = 0; column < columns; column++) { + result[offset + column] = _apply( + _data[offset + column], + vectorData[axis == Axis.columns ? column : row], + op, + ); + } + } + return Matrix32.storage(rows, columns, result, copy: false); + } + Matrix32 _broadcastVector32(Vector32 vector, ElementwiseOp op) { final vectorData = vector.unsafeValuesView; final result = Float32List(count); diff --git a/lib/src/vector.dart b/lib/src/vector.dart index 03c8238..aae7c0a 100644 --- a/lib/src/vector.dart +++ b/lib/src/vector.dart @@ -108,6 +108,14 @@ class Vector64 extends IterableBase { dynamic operator /(Object other) => _binary(other, ElementwiseOp.divide); + /// Scalar multiply that preserves the static [Vector64] type. + /// + /// `operator *` returns `dynamic` because it also dispatches to the + /// vector/matrix product; use [scale] when the result must stay + /// statically typed for chaining. + Vector64 scale(num factor) => + Vector64.storage(Kernels.scale(_data, factor.toDouble()), copy: false); + double dot(Vector64 other) { _checkLength(other); return Kernels.dot(_data, 0, other._data, 0, length); diff --git a/lib/src/vector32.dart b/lib/src/vector32.dart index 939db28..e0df39c 100644 --- a/lib/src/vector32.dart +++ b/lib/src/vector32.dart @@ -80,6 +80,14 @@ class Vector32 extends IterableBase { dynamic operator /(Object other) => _binary(other, ElementwiseOp.divide); + /// Scalar multiply that preserves the static [Vector32] type. + /// + /// `operator *` returns `dynamic` because it also dispatches to the + /// vector/matrix product; use [scale] when the result must stay + /// statically typed for chaining. + Vector32 scale(num factor) => + Vector32.storage(Kernels.scale32(_data, factor.toDouble()), copy: false); + double dot(Vector32 other) { _checkLength(other); return Kernels.dot32(_data, 0, other._data, 0, length); diff --git a/test/matrix_api_test.dart b/test/matrix_api_test.dart index bdb8649..5adbcf2 100644 --- a/test/matrix_api_test.dart +++ b/test/matrix_api_test.dart @@ -122,4 +122,67 @@ void main() { throwsArgumentError, ); }); + + test('typed arithmetic methods preserve static types', () { + final Matrix64 a = mat([ + [1, 2, 3], + [4, 5, 6], + ]); + + final Matrix64 scaled = a.scale(2); + expect(scaled.transpose.toRows(), [ + [2, 8], + [4, 10], + [6, 12], + ]); + + final Vector64 x = vec([1, 2, 3]); + final Vector64 doubled = x.scale(2); + expect(doubled.toList(), [2, 4, 6]); + + final Matrix32 a32 = mat32([ + [1, 2], + [3, 4], + ]); + expect(a32.scale(3).transpose.toRows(), [ + [3, 9], + [6, 12], + ]); + expect(vec32([1, 2]).scale(4).toList(), [4, 8]); + }); + + test('explicit axis broadcasting removes square-matrix ambiguity', () { + final a = mat([ + [1, 2], + [3, 4], + ]); + final v = vec([10, 20]); + + expect(a.addVector(v, Axis.columns).toRows(), [ + [11, 22], + [13, 24], + ]); + expect(a.addVector(v, Axis.rows).toRows(), [ + [11, 12], + [23, 24], + ]); + expect(a.subtractVector(v, Axis.columns).toRows(), [ + [-9, -18], + [-7, -16], + ]); + + final a32 = mat32([ + [1, 2], + [3, 4], + ]); + expect(a32.multiplyVector(vec32([2, 3]), Axis.rows).toRows(), [ + [2, 4], + [9, 12], + ]); + + expect( + () => a.addVector(vec([1, 2, 3]), Axis.columns), + throwsArgumentError, + ); + }); } From d81ebf8c1aebb5e928738b8070587d7b214afbb6 Mon Sep 17 00:00:00 2001 From: kelikovic Date: Sun, 13 Sep 2026 03:32:21 +0200 Subject: [PATCH 2/3] feat: CPU reverse-mode autodiff with nn layers and optimizers Add Variable (reverse-mode AD over Matrix64) with matmul/elementwise/activations/reductions, broadcast-with-reduce addRowVector, and memory-saving graph freeing in backward(). Add nn module: Module/Linear/ReLU/Tanh/Sigmoid/Sequential and SGD (with momentum)/Adam optimizers. Pure CPU on the existing SIMD kernels; no FFI. --- lib/matrices.dart | 2 + lib/src/autodiff.dart | 435 ++++++++++++++++++++++++++++++++++++++++ lib/src/nn.dart | 207 +++++++++++++++++++ test/autodiff_test.dart | 179 +++++++++++++++++ test/nn_test.dart | 91 +++++++++ 5 files changed, 914 insertions(+) create mode 100644 lib/src/autodiff.dart create mode 100644 lib/src/nn.dart create mode 100644 test/autodiff_test.dart create mode 100644 test/nn_test.dart diff --git a/lib/matrices.dart b/lib/matrices.dart index bb27ade..adf5f44 100644 --- a/lib/matrices.dart +++ b/lib/matrices.dart @@ -17,3 +17,5 @@ export 'src/iterative.dart'; export 'src/krylov.dart'; export 'src/advanced.dart'; export 'src/sparse_matrix.dart'; +export 'src/autodiff.dart'; +export 'src/nn.dart'; diff --git a/lib/src/autodiff.dart b/lib/src/autodiff.dart new file mode 100644 index 0000000..dd55342 --- /dev/null +++ b/lib/src/autodiff.dart @@ -0,0 +1,435 @@ +import 'dart:math' as math; +import 'dart:typed_data'; + +import 'matrix.dart'; +import 'types.dart'; +import 'vector.dart'; + +/// CPU reverse-mode automatic differentiation node. +/// +/// A [Variable] wraps a [Matrix64] [value] and, when it participates in +/// a differentiable expression, records a backward closure plus its +/// parents so that [backward] can walk the graph in reverse topological +/// order and accumulate gradients into every leaf's [grad]. +/// +/// The engine is pure Dart on top of the existing SIMD [Matrix64] +/// kernels — no FFI, no GPU. Gradient math is expressed with the same +/// matrix ops as the forward pass. Scalars are represented as `1x1` +/// matrices so a single code path covers losses and tensors. +/// +/// Values are treated as immutable: every op allocates a fresh +/// [Matrix64] rather than mutating storage in place, so the recorded +/// graph never observes a value changing under it. Do not mutate a +/// [value] (e.g. via `set`) after it has been used in an expression. +class Variable { + Matrix64 value; + Matrix64? grad; + bool requiresGrad; + + void Function()? _backward; + List _parents; + + Variable(this.value, {this.requiresGrad = false}) : _parents = const []; + + factory Variable.fromRows( + Iterable> rows, { + bool requiresGrad = false, + }) => Variable(Matrix64.fromRows(rows), requiresGrad: requiresGrad); + + factory Variable.scalar(num value, {bool requiresGrad = false}) => + Variable(Matrix64.full(1, 1, value), requiresGrad: requiresGrad); + + factory Variable.filled( + int rows, + int columns, + num value, { + bool requiresGrad = false, + }) => + Variable(Matrix64.full(rows, columns, value), requiresGrad: requiresGrad); + + factory Variable.zeros(int rows, int columns, {bool requiresGrad = false}) => + Variable(Matrix64.zeros(rows, columns), requiresGrad: requiresGrad); + + factory Variable.random( + int rows, + int columns, { + int? seed, + bool requiresGrad = false, + }) => Variable( + Matrix64.random(rows, columns, seed: seed), + requiresGrad: requiresGrad, + ); + + int get rows => value.rows; + + int get columns => value.columns; + + bool get isScalar => value.rows == 1 && value.columns == 1; + + /// The single element of a `1x1` variable (e.g. a loss). + double get item { + if (!isScalar) { + throw StateError( + 'item is only defined for 1x1 variables, got ($rows, $columns).', + ); + } + return value(0, 0); + } + + // --------------------------------------------------------------------- + // Autograd core. + // --------------------------------------------------------------------- + + /// Runs reverse-mode differentiation from this node. The root is + /// seeded with [seed], or with ones when this is a scalar. Gradients + /// accumulate into each participating node's [grad]; call [zeroGrad] + /// between optimization steps. + /// + /// When [retainGraph] is false (the default) every intermediate + /// (non-leaf) node's forward value and gradient are released as its + /// backward closure finishes, and the graph links are cleared. This + /// keeps peak memory flat across training steps — only leaves + /// (parameters and inputs) and this root's value survive. Pass + /// `retainGraph: true` to inspect intermediates or to call + /// [backward] again on the same graph. + void backward({Matrix64? seed, bool retainGraph = false}) { + if (!requiresGrad) { + throw StateError( + 'backward() called on a Variable with requiresGrad = false.', + ); + } + if (seed != null) { + if (seed.rows != rows || seed.columns != columns) { + throw ArgumentError( + 'Seed shape (${seed.rows}, ${seed.columns}) does not match ' + '($rows, $columns).', + ); + } + grad = seed; + } else if (isScalar) { + grad = Matrix64.full(1, 1, 1.0); + } else { + throw StateError( + 'backward() on a non-scalar Variable requires an explicit seed gradient.', + ); + } + + final ordered = []; + final visited = {}; + void visit(Variable node) { + if (!visited.add(node)) { + return; + } + for (final parent in node._parents) { + visit(parent); + } + ordered.add(node); + } + + visit(this); + for (var i = ordered.length - 1; i >= 0; i--) { + final node = ordered[i]; + node._backward?.call(); + if (!retainGraph && node._backward != null) { + // Op output: no backward closure reads a node's own forward + // value, and all of its consumers ran earlier in reverse + // order, so the buffers are safe to drop now. + if (!identical(node, this)) { + node.value = _released; + node.grad = null; + } + node._parents = const []; + node._backward = null; + } + } + } + + /// Clears the accumulated gradient on this node. + void zeroGrad() => grad = null; + + /// Returns a detached copy that does not participate in autograd. + Variable detach() => Variable(value.copy()); + + static int _noGradDepth = 0; + + /// Runs [body] with autograd disabled so no graph is recorded. + static T noGrad(T Function() body) { + _noGradDepth++; + try { + return body(); + } finally { + _noGradDepth--; + } + } + + void _setBackward(List parents, void Function() backward) { + if (_noGradDepth > 0) { + return; + } + _parents = parents; + _backward = backward; + requiresGrad = true; + } + + void _accumulate(Matrix64 contribution) { + final current = grad; + grad = current == null ? contribution : current + contribution; + } + + /// Empty matrix that replaces a freed intermediate's forward buffer. + static final Matrix64 _released = Matrix64.zeros(0, 0); + + // --------------------------------------------------------------------- + // Differentiable ops. + // --------------------------------------------------------------------- + + /// Matrix product. For `C = A @ B`: `dA = dC @ Bᵀ`, `dB = Aᵀ @ dC`. + Variable matmul(Variable other) { + final out = Variable(value.matmul(other.value)); + if (requiresGrad || other.requiresGrad) { + final a = this; + final b = other; + out._setBackward([a, b], () { + final g = out.grad!; + if (a.requiresGrad) { + a._accumulate(g.matmul(b.value.transpose)); + } + if (b.requiresGrad) { + b._accumulate(a.value.transpose.matmul(g)); + } + }); + } + return out; + } + + Variable operator +(Variable other) { + _checkSameShape(other); + final out = Variable(value + other.value); + if (requiresGrad || other.requiresGrad) { + final a = this; + final b = other; + out._setBackward([a, b], () { + final g = out.grad!; + if (a.requiresGrad) { + a._accumulate(g); + } + if (b.requiresGrad) { + b._accumulate(g); + } + }); + } + return out; + } + + Variable operator -(Variable other) { + _checkSameShape(other); + final out = Variable(value - other.value); + if (requiresGrad || other.requiresGrad) { + final a = this; + final b = other; + out._setBackward([a, b], () { + final g = out.grad!; + if (a.requiresGrad) { + a._accumulate(g); + } + if (b.requiresGrad) { + b._accumulate(g.scale(-1)); + } + }); + } + return out; + } + + /// Elementwise (Hadamard) product. `operator *` is elementwise; + /// use [matmul] for the matrix product. + Variable operator *(Variable other) { + _checkSameShape(other); + final out = Variable(value.hadamard(other.value)); + if (requiresGrad || other.requiresGrad) { + final a = this; + final b = other; + out._setBackward([a, b], () { + final g = out.grad!; + if (a.requiresGrad) { + a._accumulate(g.hadamard(b.value)); + } + if (b.requiresGrad) { + b._accumulate(g.hadamard(a.value)); + } + }); + } + return out; + } + + /// Scalar multiply. + Variable scale(num factor) { + final s = factor.toDouble(); + final out = Variable(value.scale(s)); + if (requiresGrad) { + final a = this; + out._setBackward([a], () => a._accumulate(out.grad!.scale(s))); + } + return out; + } + + Variable get transpose { + final out = Variable(value.transpose); + if (requiresGrad) { + final a = this; + out._setBackward([a], () => a._accumulate(out.grad!.transpose)); + } + return out; + } + + Variable relu() { + final out = Variable(value.mapValues((v) => v > 0 ? v : 0.0)); + if (requiresGrad) { + final a = this; + out._setBackward([a], () { + final mask = a.value.mapValues((v) => v > 0 ? 1.0 : 0.0); + a._accumulate(out.grad!.hadamard(mask)); + }); + } + return out; + } + + Variable sigmoid() { + final s = value.mapValues((v) => 1.0 / (1.0 + math.exp(-v))); + final out = Variable(s); + if (requiresGrad) { + final a = this; + out._setBackward([a], () { + final derivative = s.hadamard(s.mapValues((v) => 1.0 - v)); + a._accumulate(out.grad!.hadamard(derivative)); + }); + } + return out; + } + + Variable tanh() { + final t = value.mapValues(_tanh); + final out = Variable(t); + if (requiresGrad) { + final a = this; + out._setBackward([a], () { + final derivative = t.mapValues((v) => 1.0 - v * v); + a._accumulate(out.grad!.hadamard(derivative)); + }); + } + return out; + } + + Variable exp() { + final e = value.exp(); + final out = Variable(e); + if (requiresGrad) { + final a = this; + out._setBackward([a], () => a._accumulate(out.grad!.hadamard(e))); + } + return out; + } + + Variable log() { + final out = Variable(value.mapValues(math.log)); + if (requiresGrad) { + final a = this; + out._setBackward([a], () { + final inverse = a.value.mapValues((v) => 1.0 / v); + a._accumulate(out.grad!.hadamard(inverse)); + }); + } + return out; + } + + /// Sum of all elements, returning a `1x1` variable. + Variable sum() { + final out = Variable(Matrix64.full(1, 1, value.sum)); + if (requiresGrad) { + final a = this; + out._setBackward([a], () { + a._accumulate(Matrix64.full(a.rows, a.columns, out.grad!(0, 0))); + }); + } + return out; + } + + /// Mean of all elements, returning a `1x1` variable. + Variable mean() { + final n = value.count; + final out = Variable(Matrix64.full(1, 1, value.sum / n)); + if (requiresGrad) { + final a = this; + out._setBackward([a], () { + a._accumulate(Matrix64.full(a.rows, a.columns, out.grad!(0, 0) / n)); + }); + } + return out; + } + + /// Adds a `1 x columns` row [bias] to every row (broadcast across + /// rows). The bias gradient reduces by summing the incoming gradient + /// over rows — the canonical linear-layer bias rule. + Variable addRowVector(Variable bias) { + if (bias.rows != 1 || bias.columns != columns) { + throw ArgumentError( + 'Row bias must be (1, $columns), got (${bias.rows}, ${bias.columns}).', + ); + } + final biasVector = Vector64.storage( + bias.value.unsafeValuesView, + copy: false, + ); + final out = Variable(value.addVector(biasVector, Axis.columns)); + if (requiresGrad || bias.requiresGrad) { + final a = this; + final b = bias; + out._setBackward([a, b], () { + final g = out.grad!; + if (a.requiresGrad) { + a._accumulate(g); + } + if (b.requiresGrad) { + b._accumulate(_sumRows(g)); + } + }); + } + return out; + } + + void _checkSameShape(Variable other) { + if (rows != other.rows || columns != other.columns) { + throw ArgumentError( + 'Shape mismatch: ($rows, $columns) vs ' + '(${other.rows}, ${other.columns}).', + ); + } + } + + @override + String toString() => + 'Variable(($rows, $columns), requiresGrad: $requiresGrad)'; +} + +Matrix64 _sumRows(Matrix64 matrix) { + final columns = matrix.columns; + final data = Float64List(columns); + final flat = matrix.unsafeValuesView; + for (var row = 0; row < matrix.rows; row++) { + final offset = row * columns; + for (var column = 0; column < columns; column++) { + data[column] += flat[offset + column]; + } + } + return Matrix64.storage(1, columns, data, copy: false); +} + +double _tanh(double x) { + if (x > 20) { + return 1.0; + } + if (x < -20) { + return -1.0; + } + final e2 = math.exp(2 * x); + return (e2 - 1) / (e2 + 1); +} diff --git a/lib/src/nn.dart b/lib/src/nn.dart new file mode 100644 index 0000000..9d636c0 --- /dev/null +++ b/lib/src/nn.dart @@ -0,0 +1,207 @@ +import 'dart:math' as math; + +import 'autodiff.dart'; +import 'matrix.dart'; + +/// Base class for differentiable, parameterized layers built on +/// [Variable]. Instances are callable: `layer(input)`. +abstract class Module { + /// Trainable parameters exposed to an optimizer. + List get parameters; + + /// Forward pass. + Variable call(Variable input); + + /// Clears gradients on all parameters. + void zeroGrad() { + for (final parameter in parameters) { + parameter.zeroGrad(); + } + } +} + +/// Fully connected layer computing `input @ weight + bias`. +/// +/// `weight` is `(inFeatures, outFeatures)` initialized with Xavier +/// uniform sampling; `bias` is `(1, outFeatures)` initialized to zero +/// and broadcast across rows. +class Linear extends Module { + final Variable weight; + final Variable bias; + + Linear(int inFeatures, int outFeatures, {int? seed}) + : weight = _initWeight(inFeatures, outFeatures, seed), + bias = Variable.zeros(1, outFeatures, requiresGrad: true); + + @override + Variable call(Variable input) => input.matmul(weight).addRowVector(bias); + + @override + List get parameters => [weight, bias]; + + static Variable _initWeight(int inFeatures, int outFeatures, int? seed) { + final limit = math.sqrt(6.0 / (inFeatures + outFeatures)); + final base = Matrix64.random(inFeatures, outFeatures, seed: seed); + return Variable( + base.mapValues((value) => (value * 2 - 1) * limit), + requiresGrad: true, + ); + } +} + +class ReLU extends Module { + @override + Variable call(Variable input) => input.relu(); + + @override + List get parameters => const []; +} + +class Tanh extends Module { + @override + Variable call(Variable input) => input.tanh(); + + @override + List get parameters => const []; +} + +class Sigmoid extends Module { + @override + Variable call(Variable input) => input.sigmoid(); + + @override + List get parameters => const []; +} + +/// Runs its [layers] in order, threading the output of each into the +/// next. +class Sequential extends Module { + final List layers; + + Sequential(this.layers); + + @override + Variable call(Variable input) { + var output = input; + for (final layer in layers) { + output = layer(output); + } + return output; + } + + @override + List get parameters => [ + for (final layer in layers) ...layer.parameters, + ]; +} + +/// Base class for gradient-descent optimizers over a parameter list. +abstract class Optimizer { + List get parameters; + + /// Applies one update step from the current gradients. + void step(); + + /// Clears gradients on all parameters. + void zeroGrad() { + for (final parameter in parameters) { + parameter.zeroGrad(); + } + } +} + +/// Vanilla stochastic gradient descent with optional [momentum]. +/// [step] updates parameters in a [Variable.noGrad] scope so the update +/// is not recorded on the graph. +class SGD extends Optimizer { + @override + final List parameters; + final double learningRate; + final double momentum; + final List _velocity; + + SGD(this.parameters, {this.learningRate = 0.01, this.momentum = 0.0}) + : _velocity = List.filled(parameters.length, null); + + @override + void step() { + Variable.noGrad(() { + for (var i = 0; i < parameters.length; i++) { + final gradient = parameters[i].grad; + if (gradient == null) { + continue; + } + if (momentum > 0) { + final previous = _velocity[i]; + final velocity = previous == null + ? gradient + : previous.scale(momentum) + gradient; + _velocity[i] = velocity; + parameters[i].value = + parameters[i].value - velocity.scale(learningRate); + } else { + parameters[i].value = + parameters[i].value - gradient.scale(learningRate); + } + } + }); + } +} + +/// Adam optimizer with bias-corrected first and second moment estimates. +class Adam extends Optimizer { + @override + final List parameters; + final double learningRate; + final double beta1; + final double beta2; + final double epsilon; + + final List _firstMoment; + final List _secondMoment; + int _timestep = 0; + + Adam( + this.parameters, { + this.learningRate = 0.001, + this.beta1 = 0.9, + this.beta2 = 0.999, + this.epsilon = 1e-8, + }) : _firstMoment = List.filled(parameters.length, null), + _secondMoment = List.filled(parameters.length, null); + + @override + void step() { + _timestep++; + final firstCorrection = 1 - math.pow(beta1, _timestep).toDouble(); + final secondCorrection = 1 - math.pow(beta2, _timestep).toDouble(); + Variable.noGrad(() { + for (var i = 0; i < parameters.length; i++) { + final gradient = parameters[i].grad; + if (gradient == null) { + continue; + } + final previousFirst = _firstMoment[i]; + final first = previousFirst == null + ? gradient.scale(1 - beta1) + : previousFirst.scale(beta1) + gradient.scale(1 - beta1); + + final gradientSquared = gradient.hadamard(gradient); + final previousSecond = _secondMoment[i]; + final second = previousSecond == null + ? gradientSquared.scale(1 - beta2) + : previousSecond.scale(beta2) + gradientSquared.scale(1 - beta2); + + _firstMoment[i] = first; + _secondMoment[i] = second; + + final firstHat = first.scale(1 / firstCorrection); + final secondHat = second.scale(1 / secondCorrection); + final update = firstHat.hadamard( + secondHat.mapValues((value) => 1.0 / (math.sqrt(value) + epsilon)), + ); + parameters[i].value = parameters[i].value - update.scale(learningRate); + } + }); + } +} diff --git a/test/autodiff_test.dart b/test/autodiff_test.dart new file mode 100644 index 0000000..31b123f --- /dev/null +++ b/test/autodiff_test.dart @@ -0,0 +1,179 @@ +import 'dart:math' as math; + +import 'package:matrices/matrices.dart'; +import 'package:test/test.dart'; + +void main() { + test('matmul backward matches finite differences', () { + final a = Variable.fromRows([ + [1.0, 2.0, 3.0], + [4.0, 5.0, 6.0], + ], requiresGrad: true); + final b = Variable.fromRows([ + [0.5, -1.0], + [2.0, 1.5], + [-0.5, 0.25], + ], requiresGrad: true); + + double loss(Matrix64 av, Matrix64 bv) => av.matmul(bv).sum; + + (a.matmul(b)).sum().backward(); + + _checkGradient(a, b, loss); + }); + + test('nonlinear scalar graph matches finite differences', () { + final x = Variable.fromRows([ + [0.3, -0.7], + [1.2, 0.4], + ], requiresGrad: true); + + Variable forward(Variable v) => + v.sigmoid().matmul(v.transpose).relu().sum(); + + forward(x).backward(); + + final analytic = x.grad!; + final numeric = _numericalGradient( + x.value, + (m) => forward(Variable(m)).item, + ); + _expectClose(analytic, numeric, 1e-5); + }); + + test('linear regression converges via gradient descent', () { + final rng = math.Random(7); + const samples = 40; + const features = 3; + + final trueWeights = Variable.fromRows([ + [2.0], + [-3.0], + [0.5], + ]); + + final xRows = List>.generate( + samples, + (_) => List.generate(features, (_) => rng.nextDouble() * 2 - 1), + ); + final x = Variable.fromRows(xRows); + final y = x.matmul(trueWeights); + + final w = Variable.zeros(features, 1, requiresGrad: true); + const learningRate = 0.5; + var lastLoss = double.infinity; + + for (var step = 0; step < 400; step++) { + w.zeroGrad(); + final prediction = x.matmul(w); + final error = prediction - y; + final loss = (error * error).mean(); + loss.backward(); + + final gradient = w.grad!; + Variable.noGrad(() { + w.value = w.value - gradient.scale(learningRate); + }); + lastLoss = loss.item; + } + + expect(lastLoss, lessThan(1e-8)); + for (var i = 0; i < features; i++) { + expect(w.value(i, 0), closeTo(trueWeights.value(i, 0), 1e-3)); + } + }); + + test('noGrad suppresses graph construction', () { + final x = Variable.scalar(2.0, requiresGrad: true); + final y = Variable.noGrad(() => x.scale(3)); + expect(y.requiresGrad, isFalse); + }); + + test('addRowVector backward reduces bias gradient over rows', () { + final x = Variable.fromRows([ + [1.0, 2.0, 3.0], + [4.0, 5.0, 6.0], + ], requiresGrad: true); + final bias = Variable.fromRows([ + [0.5, -1.0, 2.0], + ], requiresGrad: true); + + x.addRowVector(bias).sum().backward(); + + // d/dx sum(x + bias) = ones like x. + _expectClose(x.grad!, Matrix64.full(2, 3, 1.0), 1e-9); + // d/dbias = column sums of the ones gradient = number of rows. + _expectClose(bias.grad!, Matrix64.full(1, 3, 2.0), 1e-9); + }); + + test('backward frees intermediate buffers by default', () { + final w = Variable.fromRows([ + [1.0], + [2.0], + ], requiresGrad: true); + final x = Variable.fromRows([ + [3.0, 4.0], + ]); + + final hidden = x.matmul(w); + final loss = hidden.sum(); + loss.backward(); + + // Intermediate value + grad released; parameter grad retained. + expect(hidden.value.count, 0); + expect(hidden.grad, isNull); + expect(w.grad, isNotNull); + }); + + test('retainGraph keeps intermediates alive', () { + final w = Variable.fromRows([ + [1.0], + [2.0], + ], requiresGrad: true); + final x = Variable.fromRows([ + [3.0, 4.0], + ]); + + final hidden = x.matmul(w); + hidden.sum().backward(retainGraph: true); + + expect(hidden.value.count, greaterThan(0)); + }); +} + +void _checkGradient( + Variable a, + Variable b, + double Function(Matrix64, Matrix64) loss, +) { + final numericA = _numericalGradient(a.value, (m) => loss(m, b.value)); + final numericB = _numericalGradient(b.value, (m) => loss(a.value, m)); + _expectClose(a.grad!, numericA, 1e-5); + _expectClose(b.grad!, numericB, 1e-5); +} + +Matrix64 _numericalGradient(Matrix64 input, double Function(Matrix64) loss) { + const epsilon = 1e-6; + final base = input.values; + final grad = List.filled(base.length, 0); + for (var i = 0; i < base.length; i++) { + final plus = base.sublist(0); + final minus = base.sublist(0); + plus[i] += epsilon; + minus[i] -= epsilon; + final lossPlus = loss(Matrix64.fromFlat(plus, input.rows, input.columns)); + final lossMinus = loss(Matrix64.fromFlat(minus, input.rows, input.columns)); + grad[i] = (lossPlus - lossMinus) / (2 * epsilon); + } + return Matrix64.fromFlat(grad, input.rows, input.columns); +} + +void _expectClose(Matrix64 actual, Matrix64 expected, double tolerance) { + expect(actual.rows, expected.rows); + expect(actual.columns, expected.columns); + for (var r = 0; r < actual.rows; r++) { + for (var c = 0; c < actual.columns; c++) { + expect(actual(r, c), closeTo(expected(r, c), tolerance)); + } + } +} diff --git a/test/nn_test.dart b/test/nn_test.dart new file mode 100644 index 0000000..c5eeb27 --- /dev/null +++ b/test/nn_test.dart @@ -0,0 +1,91 @@ +import 'package:matrices/matrices.dart'; +import 'package:test/test.dart'; + +void main() { + test('Linear forward applies weight and bias', () { + final layer = Linear(2, 3, seed: 1); + final input = Variable.fromRows([ + [1.0, 2.0], + [3.0, 4.0], + ]); + + final output = layer(input); + expect(output.rows, 2); + expect(output.columns, 3); + expect(layer.parameters, hasLength(2)); + }); + + test('MLP learns a nonlinear XOR-like mapping', () { + final model = Sequential([ + Linear(2, 8, seed: 3), + Tanh(), + Linear(8, 1, seed: 4), + ]); + final optimizer = SGD(model.parameters, learningRate: 0.1); + + final inputs = Variable.fromRows([ + [0.0, 0.0], + [0.0, 1.0], + [1.0, 0.0], + [1.0, 1.0], + ]); + final targets = Variable.fromRows([ + [0.0], + [1.0], + [1.0], + [0.0], + ]); + + var lastLoss = double.infinity; + for (var epoch = 0; epoch < 2000; epoch++) { + optimizer.zeroGrad(); + final prediction = model(inputs); + final error = prediction - targets; + final loss = (error * error).mean(); + loss.backward(); + optimizer.step(); + lastLoss = loss.item; + } + + expect(lastLoss, lessThan(1e-3)); + + final prediction = Variable.noGrad(() => model(inputs)); + expect(prediction.value(0, 0), closeTo(0, 0.1)); + expect(prediction.value(1, 0), closeTo(1, 0.1)); + expect(prediction.value(2, 0), closeTo(1, 0.1)); + expect(prediction.value(3, 0), closeTo(0, 0.1)); + }); + + test('SGD with momentum minimizes a quadratic', () { + final w = Variable.fromRows([ + [5.0, -3.0], + ], requiresGrad: true); + final optimizer = SGD([w], learningRate: 0.01, momentum: 0.9); + + for (var step = 0; step < 300; step++) { + optimizer.zeroGrad(); + // Minimize sum(w^2); minimum at the origin. + (w * w).sum().backward(); + optimizer.step(); + } + + expect(w.value(0, 0), closeTo(0, 1e-3)); + expect(w.value(0, 1), closeTo(0, 1e-3)); + }); + + test('Adam minimizes a quadratic toward the origin', () { + final w = Variable.fromRows([ + [5.0, -3.0], + ], requiresGrad: true); + final optimizer = Adam([w], learningRate: 0.1); + + for (var step = 0; step < 400; step++) { + optimizer.zeroGrad(); + (w * w).sum().backward(); + optimizer.step(); + } + + expect(w.value(0, 0), closeTo(0, 1e-3)); + expect(w.value(0, 1), closeTo(0, 1e-3)); + }); +} From 5acd17278556751671381d0d333f620115b8412f Mon Sep 17 00:00:00 2001 From: kelikovic Date: Sun, 13 Sep 2026 03:36:24 +0200 Subject: [PATCH 3/3] docs: changelog for typed ops and CPU autodiff/nn --- CHANGELOG.md | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 8e876ca..30aa6e6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,16 @@ +## Unreleased + +- Add statically-typed `scale()` to `Matrix64`, `Matrix32`, `Vector64`, and + `Vector32`, and explicit axis broadcasting (`addVector`, `subtractVector`, + `multiplyVector`, `divideVector`) on the matrix types to give a fully typed + multiply path and remove the square-matrix broadcast ambiguity. +- Add a CPU reverse-mode automatic differentiation layer (`Variable`) built on + the existing SIMD kernels, covering matmul, elementwise arithmetic, + activations, reductions, broadcast-with-reduce bias, and memory-saving graph + freeing in `backward`. +- Add a neural-network module: `Module`, `Linear`, `ReLU`, `Tanh`, `Sigmoid`, + `Sequential`, and `SGD` (with momentum) / `Adam` optimizers. + ## 2.0.0 - Redesign the package around high performance and script-friendly ergonomics,