diff --git a/AGENTS.md b/AGENTS.md
index 6ed1bf4..52a46ff 100644
--- a/AGENTS.md
+++ b/AGENTS.md
@@ -165,6 +165,38 @@ val t2dNested = Tensor2(Axis[A], Axis[B]).fromArray(
 )
 ```
 
+### Identity Matrices with `eye`
+
+`eye` is a creation method on the rank 2 factory: fix the shape first, then ask for the identity matrix.
+
+```scala
+// Identity matrix, both axes labelled explicitly
+val eye = Tensor2(Axis[A] -> 3, Axis[B] -> 3).eye
+
+// From a single extent: the second axis is the primed copy of the first,
+// i.e. the type is Tensor2[A, Prime[A], Float32]
+val primedEye = Tensor2(Axis[A] -> 3).eye
+
+// From a shape
+val eyeFromShape = Tensor2(Shape2(Axis[A] -> 3, Axis[B] -> 3)).eye
+
+// Non-square: the diagonal stops at the shorter axis
+val wideEye = Tensor2(Axis[A] -> 2, Axis[B] -> 3).eye
+
+// Unlike fill and fromArray, eye has no values to derive the value type from.
+// It defaults to Float32 and takes the value type as an argument.
+val intEye = Tensor2(Axis[A] -> 3, Axis[B] -> 3).eye(VType[Int32])
+```
+
+```scala
+// ERROR: eye only exists on the rank 2 factory
+val notAMatrix = Tensor1(Axis[A] -> 3).eye
+// error:
+// value eye is not a member of dimwit.tensor.Tensor.ShapedFactory[Tuple1[repl.MdocSession.MdocApp.A]]
+// val notAMatrix = Tensor1(Axis[A] -> 3).eye
+//                  ^^^^^^^^^^^^^^^^^^^^^^^^^
+```
+
 ### Type Aliases for Common Shapes
 
 ```scala
@@ -358,10 +390,10 @@ val wrong = t.sum(Axis[C])
 // Conflicting definitions:
 // val t:
 //   dimwit.tensor.Tensor2[MdocApp0.this.A, MdocApp0.this.B,
-//     dimwit.tensor.DType.Float32] in class MdocApp0 at line 53 and
+//     dimwit.tensor.DType.Float32] in class MdocApp0 at line 58 and
 // val t:
 //   dimwit.tensor.Tensor2[MdocApp0.this.A, MdocApp0.this.B,
-//     dimwit.tensor.DType.Float32] in class MdocApp0 at line 99
+//     dimwit.tensor.DType.Float32] in class MdocApp0 at line 104
 //
 ```
 
@@ -403,10 +435,10 @@ val wrong = t + 5.0f  // Use +! instead
 // Conflicting definitions:
 // val t:
 //   dimwit.tensor.Tensor2[MdocApp0.this.A, MdocApp0.this.B,
-//     dimwit.tensor.DType.Float32] in class MdocApp0 at line 53 and
+//     dimwit.tensor.DType.Float32] in class MdocApp0 at line 58 and
 // val t:
 //   dimwit.tensor.Tensor2[MdocApp0.this.A, MdocApp0.this.B,
-//     dimwit.tensor.DType.Float32] in class MdocApp0 at line 108
+//     dimwit.tensor.DType.Float32] in class MdocApp0 at line 113
 //
 ```
 
@@ -496,19 +528,19 @@ val wrong = m1.dot(Axis[B])(m2)
 // Conflicting definitions:
 // val m1:
 //   dimwit.tensor.Tensor2[MdocApp1.this.A, MdocApp1.this.B,
-//     dimwit.tensor.DType.Float32] in class MdocApp1 at line 130 and
+//     dimwit.tensor.DType.Float32] in class MdocApp1 at line 135 and
 // val m1:
 //   dimwit.tensor.Tensor2[MdocApp1.this.A, MdocApp1.this.B,
-//     dimwit.tensor.DType.Float32] in class MdocApp1 at line 133
+//     dimwit.tensor.DType.Float32] in class MdocApp1 at line 138
 // 
 // error: 
 // Conflicting definitions:
 // val m2:
 //   dimwit.tensor.Tensor2[MdocApp1.this.B, MdocApp1.this.C,
-//     dimwit.tensor.DType.Float32] in class MdocApp1 at line 131 and
+//     dimwit.tensor.DType.Float32] in class MdocApp1 at line 136 and
 // val m2:
 //   dimwit.tensor.Tensor2[MdocApp1.this.C, MdocApp1.this.D,
-//     dimwit.tensor.DType.Float32] in class MdocApp1 at line 134
+//     dimwit.tensor.DType.Float32] in class MdocApp1 at line 139
 //
 ```
 
@@ -1283,7 +1315,7 @@ val wrong = t.sum(Axis[C])  // Axis[C] not in tensor
 val t = Tensor2(Axis[A], Axis[B]).fill(1.0f)
 val wrong = t.vmap(Axis[C])(_.sum)  // Axis[C] doesn't exist
 // error: 
-// value fill is not a member of dimwit.tensor.Tensor2.DefaultsFactory[MdocApp12.this.A, MdocApp12.this.B]
+// value fill is not a member of dimwit.tensor.Tensor2.Axes2Factory[MdocApp12.this.A, MdocApp12.this.B]
 ```
 
 ### Gradient Errors
diff --git a/core/src/main/scala/dimwit/tensor/Tensor.scala b/core/src/main/scala/dimwit/tensor/Tensor.scala
index 027dba5..11db3d7 100644
--- a/core/src/main/scala/dimwit/tensor/Tensor.scala
+++ b/core/src/main/scala/dimwit/tensor/Tensor.scala
@@ -81,11 +81,11 @@ object Tensor:
 
   type IndicesOf[T <: Tuple] = Tuple.Map[T, [_] =>> Int]
 
-  /** Factory for createing tensors with a specific shape and default value type.
+  /** Factory for createing tensors with a specific shape.
     *
     * @param shape The shape of the tensor to create.
     */
-  case class DefaultsFactory[T <: Tuple: Labels](shape: Shape[T]):
+  class ShapedFactory[T <: Tuple: Labels](shape: Shape[T]):
 
     /** Ceates a tensor filled with the specified value. */
     def fill(value: Float): Tensor[T, Float32] = Tensor(shape, VType[Float32]).fill(value)
@@ -138,9 +138,9 @@ object Tensor:
     * @param shape The shape of the tensor to create
     * @param vtype The value type of the tensor to create
     */
-  case class TypedFactory[T <: Tuple: Labels, V](shape: Shape[T], vtype: VType[V]):
+  class ShapedTypedFactory[T <: Tuple: Labels, V](shape: Shape[T], vtype: VType[V]):
 
-    /** @see [[DefaultsFactory.fill]] */
+    /** @see [[ShapedFactory.fill]] */
     def fill(value: Boolean)(using IsBoolean[V]): Tensor[T, V] = Tensor(Jax.jnp.full(shape.dimensions.toPythonProxy, value, dtype = vtype.dtype.jaxType))
     def fill(value: Byte)(using IsInteger[V]): Tensor[T, V] = Tensor(Jax.jnp.full(shape.dimensions.toPythonProxy, value, dtype = vtype.dtype.jaxType))
     def fill(value: Short)(using IsInteger[V]): Tensor[T, V] = Tensor(Jax.jnp.full(shape.dimensions.toPythonProxy, value.toInt, dtype = vtype.dtype.jaxType))
@@ -149,7 +149,7 @@ object Tensor:
     def fill(value: Float)(using IsFloating[V]): Tensor[T, V] = Tensor(Jax.jnp.full(shape.dimensions.toPythonProxy, value, dtype = vtype.dtype.jaxType))
     def fill(value: Double)(using IsFloating[V]): Tensor[T, V] = Tensor(Jax.jnp.full(shape.dimensions.toPythonProxy, value, dtype = vtype.dtype.jaxType))
 
-    /** @see [[DefaultsFactory.fromArray]] */
+    /** @see [[ShapedFactory.fromArray]] */
     def fromArray(values: Array[Boolean])(using IsBoolean[V]): Tensor[T, V] = ArrayWriter.fromArray[T, V](shape, values)
     def fromArray(values: Array[Byte])(using IsInteger[V]): Tensor[T, V] = ArrayWriter.fromArray[T, V](shape, values)
     def fromArray(values: Array[Short])(using IsInteger[V]): Tensor[T, V] = ArrayWriter.fromArray[T, V](shape, values)
@@ -162,9 +162,9 @@ object Tensor:
     *
     * @param other The tensor to use as a template for the new tensor.
     */
-  case class LikeFactory[T <: Tuple: Labels, V](val other: Tensor[T, V]):
+  class LikeFactory[T <: Tuple: Labels, V](val other: Tensor[T, V]):
 
-    /** @see [[DefaultsFactory.fill]] */
+    /** @see [[ShapedFactory.fill]] */
     def fill(value: Boolean): Tensor[T, V] = Tensor(Jax.jnp.full(other.shape.dimensions.toPythonProxy, value, dtype = other.dtype.jaxType))
     def fill(value: Byte): Tensor[T, V] = Tensor(Jax.jnp.full(other.shape.dimensions.toPythonProxy, value, dtype = other.dtype.jaxType))
     def fill(value: Short): Tensor[T, V] = Tensor(Jax.jnp.full(other.shape.dimensions.toPythonProxy, value.toInt, dtype = other.dtype.jaxType))
@@ -173,7 +173,7 @@ object Tensor:
     def fill(value: Float): Tensor[T, V] = Tensor(Jax.jnp.full(other.shape.dimensions.toPythonProxy, value, dtype = other.dtype.jaxType))
     def fill(value: Double): Tensor[T, V] = Tensor(Jax.jnp.full(other.shape.dimensions.toPythonProxy, value, dtype = other.dtype.jaxType))
 
-    /** @see [[DefaultsFactory.fromArray]] */
+    /** @see [[ShapedFactory.fromArray]] */
     def fromArray(values: Array[Boolean])(using IsBoolean[V]): Tensor[T, V] = ArrayWriter.fromArray[T, V](other.shape, values)
     def fromArray(values: Array[Byte])(using IsInteger[V]): Tensor[T, V] = ArrayWriter.fromArray[T, V](other.shape, values)
     def fromArray(values: Array[Short])(using IsInteger[V]): Tensor[T, V] = ArrayWriter.fromArray[T, V](other.shape, values)
@@ -182,7 +182,7 @@ object Tensor:
     def fromArray(values: Array[Float])(using IsFloating[V]): Tensor[T, V] = ArrayWriter.fromArray[T, V](other.shape, values)
     def fromArray(values: Array[Double])(using IsFloating[V]): Tensor[T, V] = ArrayWriter.fromArray[T, V](other.shape, values)
 
-  /** Computes per-element values via eager Scala-side iteration. Used by [[DefaultsFactory.fromFunction]]. */
+  /** Computes per-element values via eager Scala-side iteration. Used by [[ShapedFactory.fromFunction]]. */
   private[tensor] def tabulate[T <: Tuple, V: scala.reflect.ClassTag](dims: List[Int], f: TypedIndex[T] => V): Array[V] =
     val strides = dims.scanRight(1)(_ * _).tail
     Array.tabulate(dims.product) { flatIdx =>
@@ -191,14 +191,14 @@ object Tensor:
 
   private[dimwit] def apply[T <: Tuple: Labels, V](jaxValue: Jax.PyDynamic): Tensor[T, V] = new Tensor(jaxValue)
 
-  /** Use the [[DefaultsFactory] to create a tensor */
-  def apply[T <: Tuple: Labels](shape: Shape[T]): DefaultsFactory[T] = DefaultsFactory(shape)
+  /** Use the [[ShapedFactory]] to create a tensor */
+  def apply[T <: Tuple: Labels](shape: Shape[T]): ShapedFactory[T] = ShapedFactory(shape)
 
-  /** Use the [[TypedFactory] to create a tensor */
-  def apply[T <: Tuple: Labels, V](shape: Shape[T], vtype: VType[V]): TypedFactory[T, V] = TypedFactory(shape, vtype)
+  /** Use the [[ShapedTypedFactory]] to create a tensor */
+  def apply[T <: Tuple: Labels, V](shape: Shape[T], vtype: VType[V]): ShapedTypedFactory[T, V] = ShapedTypedFactory(shape, vtype)
 
-  /** Use the [[LikeFactory] to create a tensor */
-  def like[T <: Tuple: Labels, V](template: Tensor[T, V]): Tensor.LikeFactory[T, V] = Tensor.LikeFactory(template)
+  /** Use the [[LikeFactory]] to create a tensor */
+  def like[T <: Tuple: Labels, V](template: Tensor[T, V]): LikeFactory[T, V] = LikeFactory(template)
 
 /** Type aliases for tensors of different ranks. */
 type Tensor0[V] = Tensor[EmptyTuple, V]
@@ -236,7 +236,7 @@ object Tensor0:
   given double2FloatingTensor[V: IsFloating]: Conversion[Double, Tensor0[V]] with
     def apply(value: Double): Tensor0[V] = Tensor0(VType[V])(value)
 
-  object DefaultsFactory:
+  object Value0Factory:
 
     def apply(value: Boolean): Tensor0[Bool] = Tensor0(VType[Bool])(value)
     def apply(value: Byte): Tensor0[Int8] = Tensor0(VType[Int8])(value)
@@ -246,7 +246,7 @@ object Tensor0:
     def apply(value: Float): Tensor0[Float32] = Tensor0(VType[Float32])(value)
     def apply(value: Double): Tensor0[Float64] = Tensor0(VType[Float64])(value)
 
-  case class TypedFactory[V](vtype: VType[V]):
+  class Value0TypedFactory[V](vtype: VType[V]):
 
     def apply(value: Boolean)(using IsBoolean[V]): Tensor0[V] = Tensor(Jax.jnp.array(value, dtype = vtype.dtype.jaxType))
     def apply(value: Byte)(using IsInteger[V]): Tensor0[V] = Tensor(Jax.jnp.array(value, dtype = vtype.dtype.jaxType))
@@ -256,8 +256,8 @@ object Tensor0:
     def apply(value: Float)(using IsFloating[V]): Tensor0[V] = Tensor(Jax.jnp.array(value, dtype = vtype.dtype.jaxType))
     def apply(value: Double)(using IsFloating[V]): Tensor0[V] = Tensor(Jax.jnp.array(value, dtype = vtype.dtype.jaxType))
 
-  export DefaultsFactory.*
-  def apply[V](vtype: VType[V]): TypedFactory[V] = TypedFactory(vtype)
+  export Value0Factory.*
+  def apply[V](vtype: VType[V]): Value0TypedFactory[V] = Value0TypedFactory(vtype)
 
   def like[V: Writer](template: Tensor0[V])(value: V): Tensor0[V] = Tensor(Jax.jnp.full(Shape0.dimensions.toPythonProxy, value, dtype = template.dtype.jaxType))
   def likeDType[V, T <: Tuple](template: Tensor[T, V])(value: Float): Tensor0[V] = Tensor(Jax.jnp.full(Shape0.dimensions.toPythonProxy, value, dtype = template.dtype.jaxType))
@@ -269,7 +269,7 @@ object Tensor0:
   */
 object Tensor1:
 
-  case class DefaultsFactory[L: Label](axis: Axis[L]):
+  class AxisFactory[L: Label](axis: Axis[L]):
 
     def fromArray(values: Array[Boolean]): Tensor1[L, Bool] = Tensor1(axis, VType[Bool]).fromArray(values)
     def fromArray(values: Array[Byte]): Tensor1[L, Int8] = Tensor1(axis, VType[Int8]).fromArray(values)
@@ -279,7 +279,7 @@ object Tensor1:
     def fromArray(values: Array[Float]): Tensor1[L, Float32] = Tensor1(axis, VType[Float32]).fromArray(values)
     def fromArray(values: Array[Double]): Tensor1[L, Float64] = Tensor1(axis, VType[Float64]).fromArray(values)
 
-  case class TypedFactory[L: Label, V](axis: Axis[L], vtype: VType[V]):
+  class AxisTypedFactory[L: Label, V](axis: Axis[L], vtype: VType[V]):
 
     def fromArray(values: Array[Boolean])(using IsBoolean[V]): Tensor1[L, V] = ArrayWriter.fromArray[Tuple1[L], V](Shape1(axis -> values.length), values)
     def fromArray(values: Array[Byte])(using IsInteger[V]): Tensor1[L, V] = ArrayWriter.fromArray[Tuple1[L], V](Shape1(axis -> values.length), values)
@@ -289,11 +289,11 @@ object Tensor1:
     def fromArray(values: Array[Float])(using IsFloating[V]): Tensor1[L, V] = ArrayWriter.fromArray[Tuple1[L], V](Shape1(axis -> values.length), values)
     def fromArray(values: Array[Double])(using IsFloating[V]): Tensor1[L, V] = ArrayWriter.fromArray[Tuple1[L], V](Shape1(axis -> values.length), values)
 
-  def apply[L: Label](axis: Axis[L]): DefaultsFactory[L] = DefaultsFactory(axis)
-  def apply[L: Label, V](axis: Axis[L], vtype: VType[V]): TypedFactory[L, V] = TypedFactory(axis, vtype)
+  def apply[L: Label](axis: Axis[L]): AxisFactory[L] = AxisFactory(axis)
+  def apply[L: Label, V](axis: Axis[L], vtype: VType[V]): AxisTypedFactory[L, V] = AxisTypedFactory(axis, vtype)
 
-  def apply[L: Label](axisExtent: AxisExtent[L]): Tensor.DefaultsFactory[Tuple1[L]] = Tensor.DefaultsFactory(Shape(axisExtent))
-  def apply[L: Label, V](axisExtent: AxisExtent[L], vtype: VType[V]): Tensor.TypedFactory[Tuple1[L], V] = Tensor.TypedFactory(Shape(axisExtent), vtype)
+  def apply[L: Label](axisExtent: AxisExtent[L]): Tensor.ShapedFactory[Tuple1[L]] = Tensor.ShapedFactory(Shape(axisExtent))
+  def apply[L: Label, V](axisExtent: AxisExtent[L], vtype: VType[V]): Tensor.ShapedTypedFactory[Tuple1[L], V] = Tensor.ShapedTypedFactory(Shape(axisExtent), vtype)
 
 /* Companion object for Tensors of rank 2 (matrices).
  *  Provides factory methods for creating tensors of rank 2 with various value types.
@@ -302,7 +302,7 @@ object Tensor2:
 
   type Array2D[V] = Array[Array[V]]
 
-  case class DefaultsFactory[L1: Label, L2: Label](axis1: Axis[L1], axis2: Axis[L2]):
+  class Axes2Factory[L1: Label, L2: Label](axis1: Axis[L1], axis2: Axis[L2]):
 
     def fromArray(values: Array2D[Boolean]): Tensor2[L1, L2, Bool] = Tensor2(axis1, axis2, VType[Bool]).fromArray(values)
     def fromArray(values: Array2D[Byte]): Tensor2[L1, L2, Int8] = Tensor2(axis1, axis2, VType[Int8]).fromArray(values)
@@ -312,7 +312,7 @@ object Tensor2:
     def fromArray(values: Array2D[Float]): Tensor2[L1, L2, Float32] = Tensor2(axis1, axis2, VType[Float32]).fromArray(values)
     def fromArray(values: Array2D[Double]): Tensor2[L1, L2, Float64] = Tensor2(axis1, axis2, VType[Float64]).fromArray(values)
 
-  case class TypedFactory[L1: Label, L2: Label, V](axis1: Axis[L1], axis2: Axis[L2], vtype: VType[V]):
+  class Axes2TypedFactory[L1: Label, L2: Label, V](axis1: Axis[L1], axis2: Axis[L2], vtype: VType[V]):
 
     private def createShape[V](values: Array2D[V]): Shape2[L1, L2] = Shape2(AxisExtent(axis1, values.length), AxisExtent(axis2, values.head.length))
 
@@ -324,15 +324,29 @@ object Tensor2:
     def fromArray(values: Array2D[Float])(using IsFloating[V]): Tensor2[L1, L2, V] = Tensor(createShape(values), VType[V]).fromArray(values.flatten)
     def fromArray(values: Array2D[Double])(using IsFloating[V]): Tensor2[L1, L2, V] = Tensor(createShape(values), VType[V]).fromArray(values.flatten)
 
-  def apply[L1: Label, L2: Label](axis1: Axis[L1], axis2: Axis[L2]): DefaultsFactory[L1, L2] = DefaultsFactory(axis1, axis2)
-  def apply[L1: Label, L2: Label, V](axis1: Axis[L1], axis2: Axis[L2], vtype: VType[V]): TypedFactory[L1, L2, V] = TypedFactory(axis1, axis2, vtype)
+  /** Factory for creating matrices of a known shape, i.e. the [[Tensor.ShapedFactory]] extended by
+    * the creation methods that only make sense for a matrix.
+    *
+    * @param matrixShape The shape of the matrix to create.
+    */
+  class Shaped2Factory[L1: Label, L2: Label](shape: Shape2[L1, L2]) extends Tensor.ShapedFactory[(L1, L2)](shape):
+
+    def eye: Tensor2[L1, L2, Float32] = eye(VType[Float32])
+    def eye[V](vtype: VType[V]): Tensor2[L1, L2, V] =
+      Tensor(Jax.jnp.eye(shape.dimensions(0), shape.dimensions(1), dtype = vtype.dtype.jaxType))
+
+  def apply[L1: Label, L2: Label](axis1: Axis[L1], axis2: Axis[L2]): Axes2Factory[L1, L2] = Axes2Factory(axis1, axis2)
+  def apply[L1: Label, L2: Label, V](axis1: Axis[L1], axis2: Axis[L2], vtype: VType[V]): Axes2TypedFactory[L1, L2, V] = Axes2TypedFactory(axis1, axis2, vtype)
+
+  /** Creates the factory for a square matrix, whose second axis is the primed copy of the given one,
+    * e.g. `Tensor2(Axis[A] -> 3).eye` has type `Tensor2[A, Prime[A], Float32]`.
+    */
+  def apply[L: Label](axisExtent: AxisExtent[L]): Shaped2Factory[L, Prime[L]] = Tensor2(axisExtent, Axis[Prime[L]] -> axisExtent.size)
 
-  def apply[L1: Label, L2: Label](axisExtent1: AxisExtent[L1], axisExtent2: AxisExtent[L2]): Tensor.DefaultsFactory[Tuple2[L1, L2]] = Tensor.DefaultsFactory(Shape(axisExtent1, axisExtent2))
-  def apply[L1: Label, L2: Label, V](axisExtent1: AxisExtent[L1], axisExtent2: AxisExtent[L2], vtype: VType[V]): Tensor.TypedFactory[Tuple2[L1, L2], V] = Tensor.TypedFactory(Shape(axisExtent1, axisExtent2), vtype)
+  def apply[L1: Label, L2: Label](shape: Shape2[L1, L2]): Shaped2Factory[L1, L2] = Shaped2Factory(shape)
+  def apply[L1: Label, L2: Label](axisExtent1: AxisExtent[L1], axisExtent2: AxisExtent[L2]): Shaped2Factory[L1, L2] = Tensor2(Shape2(axisExtent1, axisExtent2))
+  def apply[L1: Label, L2: Label, V](axisExtent1: AxisExtent[L1], axisExtent2: AxisExtent[L2], vtype: VType[V]): Tensor.ShapedTypedFactory[Tuple2[L1, L2], V] = Tensor.ShapedTypedFactory(Shape(axisExtent1, axisExtent2), vtype)
 
-  private def eyeImpl[L: Label, V](dim: AxisExtent[L], vtype: VType[V]): Tensor2[L, Prime[L], V] = Tensor(Jax.jnp.eye(dim.size, dtype = vtype.dtype.jaxType))
-  def eye[L: Label](dim: AxisExtent[L]): Tensor2[L, Prime[L], Float32] = eyeImpl(dim, VType[Float32])
-  def eye[L: Label, V](dim: AxisExtent[L], vtype: VType[V]): Tensor2[L, Prime[L], V] = eyeImpl(dim, vtype)
   def diag[L: Label, V](diag: Tensor1[L, V]): Tensor2[L, Prime[L], V] = Tensor(Jax.jnp.diag(diag.jaxValue))
 
 /** Companion object for Tensors of rank 3.
@@ -342,7 +356,7 @@ object Tensor3:
 
   type Array3D[V] = Array[Array[Array[V]]]
 
-  case class DefaultsFactory[L1: Label, L2: Label, L3: Label](axis1: Axis[L1], axis2: Axis[L2], axis3: Axis[L3]):
+  class Axes3Factory[L1: Label, L2: Label, L3: Label](axis1: Axis[L1], axis2: Axis[L2], axis3: Axis[L3]):
 
     def fromArray(values: Array3D[Boolean]): Tensor3[L1, L2, L3, Bool] = Tensor3(axis1, axis2, axis3, VType[Bool]).fromArray(values)
     def fromArray(values: Array3D[Byte]): Tensor3[L1, L2, L3, Int8] = Tensor3(axis1, axis2, axis3, VType[Int8]).fromArray(values)
@@ -352,7 +366,7 @@ object Tensor3:
     def fromArray(values: Array3D[Float]): Tensor3[L1, L2, L3, Float32] = Tensor3(axis1, axis2, axis3, VType[Float32]).fromArray(values)
     def fromArray(values: Array3D[Double]): Tensor3[L1, L2, L3, Float64] = Tensor3(axis1, axis2, axis3, VType[Float64]).fromArray(values)
 
-  case class TypedFactory[L1: Label, L2: Label, L3: Label, V](axis1: Axis[L1], axis2: Axis[L2], axis3: Axis[L3], vtype: VType[V]):
+  class Axes3TypedFactory[L1: Label, L2: Label, L3: Label, V](axis1: Axis[L1], axis2: Axis[L2], axis3: Axis[L3], vtype: VType[V]):
 
     private def createShape[V](values: Array3D[V]): Shape3[L1, L2, L3] = Shape3(AxisExtent(axis1, values.length), AxisExtent(axis2, values.head.length), AxisExtent(axis3, values.head.head.length))
     def fromArray(values: Array3D[Boolean])(using IsBoolean[V]): Tensor3[L1, L2, L3, V] = Tensor(createShape(values), VType[V]).fromArray(values.flatten.flatten)
@@ -363,8 +377,8 @@ object Tensor3:
     def fromArray(values: Array3D[Float])(using IsFloating[V]): Tensor3[L1, L2, L3, V] = Tensor(createShape(values), VType[V]).fromArray(values.flatten.flatten)
     def fromArray(values: Array3D[Double])(using IsFloating[V]): Tensor3[L1, L2, L3, V] = Tensor(createShape(values), VType[V]).fromArray(values.flatten.flatten)
 
-  def apply[L1: Label, L2: Label, L3: Label](axis1: Axis[L1], axis2: Axis[L2], axis3: Axis[L3]): DefaultsFactory[L1, L2, L3] = DefaultsFactory(axis1, axis2, axis3)
-  def apply[L1: Label, L2: Label, L3: Label, V](axis1: Axis[L1], axis2: Axis[L2], axis3: Axis[L3], vtype: VType[V]): TypedFactory[L1, L2, L3, V] = TypedFactory(axis1, axis2, axis3, vtype)
+  def apply[L1: Label, L2: Label, L3: Label](axis1: Axis[L1], axis2: Axis[L2], axis3: Axis[L3]): Axes3Factory[L1, L2, L3] = Axes3Factory(axis1, axis2, axis3)
+  def apply[L1: Label, L2: Label, L3: Label, V](axis1: Axis[L1], axis2: Axis[L2], axis3: Axis[L3], vtype: VType[V]): Axes3TypedFactory[L1, L2, L3, V] = Axes3TypedFactory(axis1, axis2, axis3, vtype)
 
-  def apply[L1: Label, L2: Label, L3: Label](axisExtent1: AxisExtent[L1], axisExtent2: AxisExtent[L2], axisExtent3: AxisExtent[L3]): Tensor.DefaultsFactory[Tuple3[L1, L2, L3]] = Tensor.DefaultsFactory(Shape(axisExtent1, axisExtent2, axisExtent3))
-  def apply[L1: Label, L2: Label, L3: Label, V](axisExtent1: AxisExtent[L1], axisExtent2: AxisExtent[L2], axisExtent3: AxisExtent[L3], vtype: VType[V]): Tensor.TypedFactory[Tuple3[L1, L2, L3], V] = Tensor.TypedFactory(Shape(axisExtent1, axisExtent2, axisExtent3), vtype)
+  def apply[L1: Label, L2: Label, L3: Label](axisExtent1: AxisExtent[L1], axisExtent2: AxisExtent[L2], axisExtent3: AxisExtent[L3]): Tensor.ShapedFactory[Tuple3[L1, L2, L3]] = Tensor.ShapedFactory(Shape(axisExtent1, axisExtent2, axisExtent3))
+  def apply[L1: Label, L2: Label, L3: Label, V](axisExtent1: AxisExtent[L1], axisExtent2: AxisExtent[L2], axisExtent3: AxisExtent[L3], vtype: VType[V]): Tensor.ShapedTypedFactory[Tuple3[L1, L2, L3], V] = Tensor.ShapedTypedFactory(Shape(axisExtent1, axisExtent2, axisExtent3), vtype)
diff --git a/core/src/test/scala/dimwit/autodiff/AutodiffSuite.scala b/core/src/test/scala/dimwit/autodiff/AutodiffSuite.scala
index 3de3295..beaeebe 100644
--- a/core/src/test/scala/dimwit/autodiff/AutodiffSuite.scala
+++ b/core/src/test/scala/dimwit/autodiff/AutodiffSuite.scala
@@ -70,7 +70,7 @@ class AutodiffSuite extends DimwitTest:
         val jf = Autodiff.jacobian(f)
 
         val x = Tensor1(Axis[A]).fromArray(Array(1.0f, 1.0f))
-        jf(x) should approxEqual(Tensor2.eye(x.extent(Axis[A]), x.vtype) *! 2.0f)
+        jf(x) should approxEqual(Tensor2(x.extent(Axis[A])).eye(x.vtype) *! 2.0f)
 
   describe("jacRev"):
 
@@ -84,8 +84,8 @@ class AutodiffSuite extends DimwitTest:
 
       // the first output is x2, so it depends on x2 only, and the other way round
       x1_dx1 should approxEqual(Tensor.like(x1_dx1).fill(0f))
-      x1_dx2 should approxEqual(Tensor2.eye(x1.extent(Axis[A]), x1.vtype))
-      x2_dx1 should approxEqual(Tensor2.eye(x2.extent(Axis[A]), x2.vtype))
+      x1_dx2 should approxEqual(Tensor2(x1.extent(Axis[A])).eye(x1.vtype))
+      x2_dx1 should approxEqual(Tensor2(x2.extent(Axis[A])).eye(x2.vtype))
       x2_dx2 should approxEqual(Tensor.like(x2_dx2).fill(0f))
 
     it("d¹ of f: Tensor1[A] => Tensor1[B] keeps the output axis first"):
@@ -94,7 +94,7 @@ class AutodiffSuite extends DimwitTest:
 
       val x = Tensor1(Axis[A]).fromArray(Array(1.0f, 1.0f))
       df(x).axes shouldBe List("B", "A")
-      df(x) should approxEqual((Tensor2.eye(x.extent(Axis[A])) *! 2.0f).relabelAll((Axis[B], Axis[A])))
+      df(x) should approxEqual((Tensor2(x.extent(Axis[A])).eye *! 2.0f).relabelAll((Axis[B], Axis[A])))
 
     it("d² of f(x1, x2) = sum(x1 * x2)"):
       def f(x1: Tensor1[A, Float32], x2: Tensor1[A, Float32]): Tensor0[Float32] = (x1 * x2).sum
@@ -106,8 +106,8 @@ class AutodiffSuite extends DimwitTest:
 
       // d²/dx1² and d²/dx2² vanish, the mixed partials are the identity
       x1_dx1 should approxEqual(Tensor.like(x1_dx1).fill(0f))
-      x1_dx2 should approxEqual(Tensor2.eye(x1.extent(Axis[A]), x1.vtype))
-      x2_dx1 should approxEqual(Tensor2.eye(x2.extent(Axis[A]), x2.vtype))
+      x1_dx2 should approxEqual(Tensor2(x1.extent(Axis[A])).eye(x1.vtype))
+      x2_dx1 should approxEqual(Tensor2(x2.extent(Axis[A])).eye(x2.vtype))
       x2_dx2 should approxEqual(Tensor.like(x2_dx2).fill(0f))
 
   describe("jacFwd"):
@@ -122,8 +122,8 @@ class AutodiffSuite extends DimwitTest:
 
       // the first output is x2, so it depends on x2 only, and the other way round
       x1_dx1 should approxEqual(Tensor.like(x1_dx1).fill(0f))
-      x1_dx2 should approxEqual(Tensor2.eye(x1.extent(Axis[A]), x1.vtype))
-      x2_dx1 should approxEqual(Tensor2.eye(x2.extent(Axis[A]), x2.vtype))
+      x1_dx2 should approxEqual(Tensor2(x1.extent(Axis[A])).eye(x1.vtype))
+      x2_dx1 should approxEqual(Tensor2(x2.extent(Axis[A])).eye(x2.vtype))
       x2_dx2 should approxEqual(Tensor.like(x2_dx2).fill(0f))
 
     it("d¹ of f: Tensor1[A] => Tensor1[B] keeps the output axis first"):
@@ -132,7 +132,7 @@ class AutodiffSuite extends DimwitTest:
 
       val x = Tensor1(Axis[A]).fromArray(Array(1.0f, 1.0f))
       df(x).axes shouldBe List("B", "A")
-      df(x) should approxEqual((Tensor2.eye(x.extent(Axis[A])) *! 2.0f).relabelAll((Axis[B], Axis[A])))
+      df(x) should approxEqual((Tensor2(x.extent(Axis[A])).eye *! 2.0f).relabelAll((Axis[B], Axis[A])))
 
     it("d² of f(x1, x2) = sum(x1 * x2)"):
       def f(x1: Tensor1[A, Float32], x2: Tensor1[A, Float32]): Tensor0[Float32] = (x1 * x2).sum
@@ -144,8 +144,8 @@ class AutodiffSuite extends DimwitTest:
 
       // d²/dx1² and d²/dx2² vanish, the mixed partials are the identity
       x1_dx1 should approxEqual(Tensor.like(x1_dx1).fill(0f))
-      x1_dx2 should approxEqual(Tensor2.eye(x1.extent(Axis[A]), x1.vtype))
-      x2_dx1 should approxEqual(Tensor2.eye(x2.extent(Axis[A]), x2.vtype))
+      x1_dx2 should approxEqual(Tensor2(x1.extent(Axis[A])).eye(x1.vtype))
+      x2_dx1 should approxEqual(Tensor2(x2.extent(Axis[A])).eye(x2.vtype))
       x2_dx2 should approxEqual(Tensor.like(x2_dx2).fill(0f))
 
   describe("hessian"):
@@ -162,7 +162,7 @@ class AutodiffSuite extends DimwitTest:
         val hf = Autodiff.hessian(f)
 
         val x = Tensor1(Axis[A]).fromArray(Array(1.0f, 5.0f))
-        hf(x) should approxEqual(Tensor2.eye(x.extent(Axis[A]), x.vtype) *! 2.0f)
+        hf(x) should approxEqual(Tensor2(x.extent(Axis[A])).eye(x.vtype) *! 2.0f)
 
       it("Hessian of f(x1, x2) = sum(x1 * x2)"):
         def f(x1: Tensor1[A, Float32], x2: Tensor1[A, Float32]): Tensor0[Float32] = (x1 * x2).sum
@@ -174,8 +174,8 @@ class AutodiffSuite extends DimwitTest:
         val (x1_dx1, x1_dx2) = x1Grad
         val (x2_dx1, x2_dx2) = x2Grad
         x1_dx1 should approxEqual(Tensor.like(x1_dx1).fill(0f))
-        x1_dx2 should approxEqual(Tensor2.eye(x1.extent(Axis[A]), x1.vtype) *! Tensor0(1.0f))
-        x2_dx1 should approxEqual(Tensor2.eye(x2.extent(Axis[A]), x2.vtype) *! Tensor0(1.0f))
+        x1_dx2 should approxEqual(Tensor2(x1.extent(Axis[A])).eye(x1.vtype) *! Tensor0(1.0f))
+        x2_dx1 should approxEqual(Tensor2(x2.extent(Axis[A])).eye(x2.vtype) *! Tensor0(1.0f))
         x2_dx2 should approxEqual(Tensor.like(x2_dx2).fill(0f))
 
   describe("jacobian of a function whose input and output axes differ"):
@@ -186,7 +186,7 @@ class AutodiffSuite extends DimwitTest:
 
       val x = Tensor1(Axis[A]).fromArray(Array(1.0f, 1.0f))
       jf(x).axes shouldBe List("B", "A")
-      jf(x) should approxEqual((Tensor2.eye(x.extent(Axis[A])) *! 2.0f).relabelAll((Axis[B], Axis[A])))
+      jf(x) should approxEqual((Tensor2(x.extent(Axis[A])).eye *! 2.0f).relabelAll((Axis[B], Axis[A])))
 
     it("primes an input axis that collides with an output axis"):
       def f(x: Tensor2[A, B, Float32]): Tensor1[B, Float32] = x.sum(Axis[A])
@@ -255,19 +255,19 @@ class AutodiffSuite extends DimwitTest:
       val jf = Autodiff.jacobian(f)
       val jac = jf(params)
 
-      jac.w.w should approxEqual(Tensor2.eye(params.w.extent(Axis[A])))
+      jac.w.w should approxEqual(Tensor2(params.w.extent(Axis[A])).eye)
       jac.w.b should approxEqual(Tensor.like(jac.w.b).fill(0f))
       jac.b.w should approxEqual(Tensor.like(jac.b.w).fill(0f))
-      jac.b.b should approxEqual(Tensor2.eye(params.b.extent(Axis[B])))
+      jac.b.b should approxEqual(Tensor2(params.b.extent(Axis[B])).eye)
 
     it("takes the hessian of a scalar loss over a case class tree"):
       def loss(p: JacParams): Tensor0[Float32] = (p.w * p.w).sum + (p.b * p.b).sum
       val hf = Autodiff.hessian(loss)
       val hess = hf(params)
 
-      hess.w.w should approxEqual(Tensor2.eye(params.w.extent(Axis[A])) *! 2.0f)
+      hess.w.w should approxEqual(Tensor2(params.w.extent(Axis[A])).eye *! 2.0f)
       hess.w.b should approxEqual(Tensor.like(hess.w.b).fill(0f))
-      hess.b.b should approxEqual(Tensor2.eye(params.b.extent(Axis[B])) *! 2.0f)
+      hess.b.b should approxEqual(Tensor2(params.b.extent(Axis[B])).eye *! 2.0f)
 
     it("differentiates a function returning a named tuple"):
       def f(x: Tensor1[A, Float32]): (u: Tensor1[A, Float32], v: Tensor1[A, Float32]) =
@@ -276,8 +276,8 @@ class AutodiffSuite extends DimwitTest:
 
       val x = Tensor1(Axis[A]).fromArray(Array(1.0f, 1.0f))
       val jac = jf(x)
-      jac.u should approxEqual(Tensor2.eye(x.extent(Axis[A])) *! 2.0f)
-      jac.v should approxEqual(Tensor2.eye(x.extent(Axis[A])) *! 3.0f)
+      jac.u should approxEqual(Tensor2(x.extent(Axis[A])).eye *! 2.0f)
+      jac.v should approxEqual(Tensor2(x.extent(Axis[A])).eye *! 3.0f)
 
   describe("Complex application"):
     it("case class support"):
diff --git a/core/src/test/scala/dimwit/tensor/TensorCreationSuite.scala b/core/src/test/scala/dimwit/tensor/TensorCreationSuite.scala
index b2023e8..3a3a237 100644
--- a/core/src/test/scala/dimwit/tensor/TensorCreationSuite.scala
+++ b/core/src/test/scala/dimwit/tensor/TensorCreationSuite.scala
@@ -92,3 +92,40 @@ class TensorCreationSuite extends DimwitTest:
         idx(Axis[A]).toFloat
       }
       result shouldEqual Tensor1(Axis[A]).fromArray(Array(0.0f, 1.0f, 2.0f, 3.0f))
+
+  describe("eye"):
+
+    it("square: from two extents or from a shape"):
+      val expected = Tensor2(Axis[A], Axis[B]).fromArray(
+        Array(Array(1.0f, 0.0f), Array(0.0f, 1.0f))
+      )
+      Tensor2(Axis[A] -> 2, Axis[B] -> 2).eye shouldEqual expected
+      Tensor2(Shape2(Axis[A] -> 2, Axis[B] -> 2)).eye shouldEqual expected
+
+    it("square: from a single extent, the second axis is the primed copy of the first"):
+      val result = Tensor2(Axis[A] -> 3).eye
+      result.shape shouldEqual Shape2(Axis[A] -> 3, Axis[Prime[A]] -> 3)
+      result shouldEqual Tensor2(Axis[A], Axis[B]).fromArray(
+        Array(Array(1.0f, 0.0f, 0.0f), Array(0.0f, 1.0f, 0.0f), Array(0.0f, 0.0f, 1.0f))
+      )
+
+    it("wide: more columns than rows, zero padded"):
+      val expected = Tensor2(Axis[A], Axis[B]).fromArray(
+        Array(Array(1.0f, 0.0f, 0.0f), Array(0.0f, 1.0f, 0.0f))
+      )
+      val result = Tensor2(Axis[A] -> 2, Axis[B] -> 3).eye
+      result.shape shouldEqual Shape2(Axis[A] -> 2, Axis[B] -> 3)
+      result shouldEqual expected
+
+    it("tall: more rows than columns, truncating"):
+      val expected = Tensor2(Axis[A], Axis[B]).fromArray(
+        Array(Array(1.0f, 0.0f), Array(0.0f, 1.0f), Array(0.0f, 0.0f))
+      )
+      val result = Tensor2(Axis[A] -> 3, Axis[B] -> 2).eye
+      result.shape shouldEqual Shape2(Axis[A] -> 3, Axis[B] -> 2)
+      result shouldEqual expected
+
+    it("defaults to Float32 and takes the vtype as an argument"):
+      Tensor2(Axis[A] -> 2, Axis[B] -> 3).eye.dtype shouldBe DType.Float32
+      Tensor2(Axis[A] -> 2, Axis[B] -> 3).eye(VType[Int32]).dtype shouldBe DType.Int32
+      Tensor2(Shape2(Axis[A] -> 2, Axis[B] -> 3)).eye(VType[Int16]).dtype shouldBe DType.Int16
diff --git a/docs/quickstart.md b/docs/quickstart.md
index c1694b8..790cef1 100644
--- a/docs/quickstart.md
+++ b/docs/quickstart.md
@@ -169,7 +169,7 @@ val vector = Tensor1(Axis[Feature]).fromArray(Array(1.0f, 2.0f))
 A `Tensor2` represents a matrix. The Tensor2 factory provides convenient methods to create special matrices, such as for example the identity 
 matrix:
 ```scala
-val eye = Tensor2.eye(Axis[Feature] -> 3)
+val eye = Tensor2(Axis[Feature] -> 3, Axis[Batch] -> 3).eye
 ``` 
 Of course, we can also create a Tensor2 from an array of data, just like we did for the general Tensor factory:
 ```scala
diff --git a/mdocs/AGENTS.md b/mdocs/AGENTS.md
index f009c64..4b6e003 100644
--- a/mdocs/AGENTS.md
+++ b/mdocs/AGENTS.md
@@ -133,6 +133,34 @@ val t2dNested = Tensor2(Axis[A], Axis[B]).fromArray(
 )
 ```
 
+### Identity Matrices with `eye`
+
+`eye` is a creation method on the rank 2 factory: fix the shape first, then ask for the identity matrix.
+
+```scala mdoc:silent
+// Identity matrix, both axes labelled explicitly
+val eye = Tensor2(Axis[A] -> 3, Axis[B] -> 3).eye
+
+// From a single extent: the second axis is the primed copy of the first,
+// i.e. the type is Tensor2[A, Prime[A], Float32]
+val primedEye = Tensor2(Axis[A] -> 3).eye
+
+// From a shape
+val eyeFromShape = Tensor2(Shape2(Axis[A] -> 3, Axis[B] -> 3)).eye
+
+// Non-square: the diagonal stops at the shorter axis
+val wideEye = Tensor2(Axis[A] -> 2, Axis[B] -> 3).eye
+
+// Unlike fill and fromArray, eye has no values to derive the value type from.
+// It defaults to Float32 and takes the value type as an argument.
+val intEye = Tensor2(Axis[A] -> 3, Axis[B] -> 3).eye(VType[Int32])
+```
+
+```scala mdoc:fail
+// ERROR: eye only exists on the rank 2 factory
+val notAMatrix = Tensor1(Axis[A] -> 3).eye
+```
+
 ### Type Aliases for Common Shapes
 
 ```scala mdoc:silent
diff --git a/mdocs/docs/quickstart.md b/mdocs/docs/quickstart.md
index 05e7dbd..4037d30 100644
--- a/mdocs/docs/quickstart.md
+++ b/mdocs/docs/quickstart.md
@@ -165,7 +165,7 @@ val vector = Tensor1(Axis[Feature]).fromArray(Array(1.0f, 2.0f))
 A `Tensor2` represents a matrix. The Tensor2 factory provides convenient methods to create special matrices, such as for example the identity 
 matrix:
 ```scala mdoc:silent
-val eye = Tensor2.eye(Axis[Feature] -> 3)
+val eye = Tensor2(Axis[Feature] -> 3, Axis[Batch] -> 3).eye
 ``` 
 Of course, we can also create a Tensor2 from an array of data, just like we did for the general Tensor factory:
 ```scala mdoc:silent