diff --git a/docs/_static/matrix.css b/docs/_static/matrix.css
new file mode 100644
index 000000000..5700ea3fc
--- /dev/null
+++ b/docs/_static/matrix.css
@@ -0,0 +1,104 @@
+/* Based on the stylesheet used by matrepr (https://github.com/alugowski/matrepr) and modified for sphinx */
+
+table.matrix  {
+    border-collapse: collapse;
+    border: 0px;
+}
+
+/* Disable a horizintal line from the default stylesheet */
+.table.matrix > :not(caption) > * > * {
+    border-bottom-width: 0px;
+}
+
+/* row indices */
+table.matrix > tbody tr th {
+    font-size: smaller;
+    font-weight: bolder;
+    vertical-align: middle;
+    text-align: right;
+}
+/* row indices are often made bold in the source data; here make them match the boldness of the th column label style*/
+table.matrix strong {
+    font-weight: bold;
+}
+
+/* column indices */
+table.matrix > thead tr th {
+    font-size: smaller;
+    font-weight: bolder;
+    vertical-align: middle;
+    text-align: center;
+}
+
+/* cells */
+table.matrix > tbody tr td {
+    vertical-align: middle;
+    text-align: center;
+    position: relative;
+}
+
+/* left border */
+table.matrix > tbody tr td:first-of-type {
+    border-left: solid 2px var(--pst-color-text-base);
+}
+/* right border */
+table.matrix > tbody tr td:last-of-type {
+    border-right: solid 2px var(--pst-color-text-base);
+}
+
+/* prevents empty cells from collapsing, especially empty rows */
+table.matrix > tbody tr td:empty::before {
+    /* basicaly fills empty cells with   */
+    content: "\00a0\00a0\00a0";
+    visibility: hidden;
+}
+table.matrix > tbody tr td:empty::after {
+    content: "\00a0\00a0\00a0";
+    visibility: hidden;
+}
+
+/* matrix bracket ticks */
+table.matrix > tbody > tr:first-child > td:first-of-type::before {
+    content: "";
+    width: 4px;
+    position: absolute;
+    top: 0;
+    bottom: 0;
+    visibility: visible;
+    left: 0;
+    right: auto;
+    border-top: solid 2px var(--pst-color-text-base);
+}
+table.matrix > tbody > tr:last-child > td:first-of-type::before {
+    content: "";
+    width: 4px;
+    position: absolute;
+    top: 0;
+    bottom: 0;
+    visibility: visible;
+    left: 0;
+    right: auto;
+    border-bottom: solid 2px var(--pst-color-text-base);
+}
+table.matrix > tbody > tr:first-child > td:last-of-type::after {
+    content: "";
+    width: 4px;
+    position: absolute;
+    top: 0;
+    bottom: 0;
+    visibility: visible;
+    left: auto;
+    right: 0;
+    border-top: solid 2px var(--pst-color-text-base);
+}
+table.matrix > tbody > tr:last-child > td:last-of-type::after {
+    content: "";
+    width: 4px;
+    position: absolute;
+    top: 0;
+    bottom: 0;
+    visibility: visible;
+    left: auto;
+    right: 0;
+    border-bottom: solid 2px var(--pst-color-text-base);
+}
diff --git a/docs/conf.py b/docs/conf.py
index 07a373203..2e6f616d8 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -36,7 +36,7 @@
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
 # ones.
 extensions = ["sphinx.ext.autodoc", "numpydoc", "sphinx_panels", "nbsphinx"]
-html_css_files = ["custom.css"]
+html_css_files = ["custom.css", "matrix.css"]
 html_js_files = ["custom.js"]
 
 # Add any paths that contain templates here, relative to this directory.
diff --git a/docs/user_guide/io.rst b/docs/user_guide/io.rst
index c13fda5d6..ecb4c0862 100644
--- a/docs/user_guide/io.rst
+++ b/docs/user_guide/io.rst
@@ -29,6 +29,7 @@ array will match the collection dtype.
     v = gb.Vector.from_coo([1, 3, 6], [2, 3, 4], float, size=10)
 
 .. csv-table::
+    :class: matrix
     :header: 0,1,2,3,4,5,6,7,8,9,10
 
     ,2.0,,3.0,,,4.0,,,
diff --git a/docs/user_guide/operations.rst b/docs/user_guide/operations.rst
index ede2efb06..3f710dc23 100644
--- a/docs/user_guide/operations.rst
+++ b/docs/user_guide/operations.rst
@@ -45,8 +45,9 @@ a Vector is treated as an nx1 column matrix.
     C << gb.semiring.min_plus(A @ B)  # functional style
 
 .. csv-table:: A
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2,3
+    :stub-columns: 1
 
     **0**,,2.0,5.0,
     **1**,,,1.5,4.25
@@ -54,8 +55,9 @@ a Vector is treated as an nx1 column matrix.
     **3**,,,,
 
 .. csv-table:: B
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,,3.0,2.0
     **1**,9.0,6.0,
@@ -63,8 +65,9 @@ a Vector is treated as an nx1 column matrix.
     **3**,0.0,5.0,
 
 .. csv-table:: C << min_plus(A @ B)
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,11.0,8.0,6.0
     **1**,4.25,4.5,2.5
@@ -90,8 +93,9 @@ a Vector is treated as an nx1 column matrix.
     w << gb.semiring.plus_times(A @ v)  # functional style
 
 .. csv-table:: A
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2,3
+    :stub-columns: 1
 
     **0**,,2.0,5.0,
     **1**,,,1.5,4.25
@@ -99,13 +103,13 @@ a Vector is treated as an nx1 column matrix.
     **3**,,,,
 
 .. csv-table:: v
-    :class: inline
+    :class: inline matrix
     :header: 0,1,2,3
 
     10.0,20.0,,40.0
 
 .. csv-table:: w << plus_times(A @ v)
-    :class: inline
+    :class: inline matrix
     :header: 0,1,2,3
 
     40.0,170.0,20.0,
@@ -127,14 +131,15 @@ a Vector is treated as an nx1 column matrix.
     u << gb.semiring.plus_plus(v @ B)  # functional style
 
 .. csv-table:: v
-    :class: inline
+    :class: inline matrix
     :header: 0,1,2,3
 
     10.0,20.0,,40.0
 
 .. csv-table:: B
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,,3.0,2.0
     **1**,9.0,6.0,
@@ -142,7 +147,7 @@ a Vector is treated as an nx1 column matrix.
     **3**,0.0,5.0,
 
 .. csv-table:: u << plus_plus(v @ B)
-    :class: inline
+    :class: inline matrix
     :header: 0,1,2
 
     69.0,84.0,12.0
@@ -181,24 +186,27 @@ Example usage:
     C << gb.binary.min(A & B)  # functional style
 
 .. csv-table:: A
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,,2.0,5.0
     **1**,1.5,,4.0
     **2**,,0.5,
 
 .. csv-table:: B
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,,3.0,-2.0
     **1**,0.0,6.0,
     **2**,,3.0,1.0
 
 .. csv-table:: C << min(A & B)
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,,2.0,-2.0
     **1**,0.0,,
@@ -265,24 +273,27 @@ should be used with the functional syntax, ``left_default`` and ``right_default`
     C << gb.binary.minus(A | B)  # functional style
 
 .. csv-table:: A
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,9.0,2.0,5.0
     **1**,1.5,,4.0
     **2**,,,
 
 .. csv-table:: B
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,4.0,0.0,-2.0
     **1**,,,
     **2**,6.0,3.0,1.0
 
 .. csv-table:: C << A.ewise_add(B, 'minus')
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2,
+    :stub-columns: 1
 
     **0**,5.0,2.0,7.0
     **1**,1.5,,4.0
@@ -310,24 +321,27 @@ should be used with the functional syntax, ``left_default`` and ``right_default`
     C << gb.binary.minus(A | B, left_default=0, right_default=0)  # functional style
 
 .. csv-table:: A
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,9.0,2.0,5.0
     **1**,1.5,,4.0
     **2**,,,
 
 .. csv-table:: B
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,4.0,0.0,-2.0
     **1**,,,
     **2**,6.0,3.0,1.0
 
 .. csv-table:: C << A.ewise_union(B, 'minus', 0, 0)
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2,
+    :stub-columns: 1
 
     **0**,5.0,2.0,7.0
     **1**,1.5,,4.0
@@ -362,13 +376,13 @@ Vector Slice Example:
     w << v[:4]
 
 .. csv-table:: v
-    :class: inline
+    :class: inline matrix
     :header: 0,1,2,3,4,5,6
 
     10.0,2.0,,40.0,-5.0,,24.0
 
 .. csv-table:: w << v[:4]
-    :class: inline
+    :class: inline matrix
     :header: 0,1,2,3
 
     10.0,2.0,,40.0
@@ -387,16 +401,18 @@ Matrix List Example:
     C << A[[0, 2], :]
 
 .. csv-table:: A
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,,2.0,5.0
     **1**,1.5,4.0,
     **2**,0.5,,-7.0
 
 .. csv-table:: C << A[[0, 2], :]
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,,2.0,5.0
     **1**,0.5,,-7.0
@@ -434,23 +450,26 @@ Matrix-Matrix Assignment Example:
     A[::2, ::2] << B
 
 .. csv-table:: A
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,,2.0,5.0
     **1**,1.5,4.0,
     **2**,0.5,,-7.0
 
 .. csv-table:: B
-    :class: inline
+    :class: inline matrix
     :header: ,0,1
+    :stub-columns: 1
 
     **0**,-99.0,-98.0
     **1**,-97.0,-96.0
 
 .. csv-table:: A[::2, ::2] << B
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,-99.0,2.0,-98.0
     **1**,1.5,4.0,
@@ -470,22 +489,24 @@ Matrix-Vector Assignment Example:
     A[1, :] << v
 
 .. csv-table:: A
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,,2.0,5.0
     **1**,1.5,4.0,
     **2**,0.5,,-7.0
 
 .. csv-table:: v
-    :class: inline
+    :class: inline matrix
     :header: 0,1,2
 
     ,,-99.0
 
 .. csv-table:: A[1, :] << v
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,,2.0,5.0
     **1**,,,-99.0
@@ -500,13 +521,13 @@ Vector-Scalar Assignment Example:
     v[:4] << 99
 
 .. csv-table:: v
-    :class: inline
+    :class: inline matrix
     :header: 0,1,2,3,4,5,6
 
     10,2,,40,-5,,24
 
 .. csv-table:: v[:4] << 99
-    :class: inline
+    :class: inline matrix
     :header: 0,1,2,3,4,5,6
 
     99,99,99,99,-5,,24
@@ -535,13 +556,13 @@ function with the collection as the argument.
     w << gb.unary.minv(v)
 
 .. csv-table:: v
-    :class: inline
+    :class: inline matrix
     :header: 0,1,2,3
 
     10.0,20.0,,40.0
 
 .. csv-table:: w << minv(v)
-    :class: inline
+    :class: inline matrix
     :header: 0,1,2,3
 
     0.1,0.05,,0.025
@@ -558,13 +579,13 @@ function with the collection as the argument.
     w << gb.indexunary.index(v)
 
 .. csv-table:: v
-    :class: inline
+    :class: inline matrix
     :header: 0,1,2,3
 
     10.0,20.0,,40.0
 
 .. csv-table:: w << index(v)
-    :class: inline
+    :class: inline matrix
     :header: 0,1,2,3
 
     0,1,,3
@@ -582,13 +603,13 @@ function with the collection as the argument.
     w << v - 15
 
 .. csv-table:: v
-    :class: inline
+    :class: inline matrix
     :header: 0,1,2,3
 
     10.0,20.0,,40.0
 
 .. csv-table:: w << v.apply('minus', right=15)
-    :class: inline
+    :class: inline matrix
     :header: 0,1,2,3,
 
     -5.0,5.0,,25.0
@@ -616,16 +637,18 @@ Upper Triangle Example:
     C << gb.select.triu(A)
 
 .. csv-table:: A
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,,2.0,5.0
     **1**,1.5,,4.0
     **2**,,0.5,-7.0
 
 .. csv-table:: C << select.triu(A)
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,,2.0,5.0
     **1**,,,4.0
@@ -643,13 +666,13 @@ Select by Value Example:
     w << gb.select.value(v >= 5)
 
 .. csv-table:: v
-    :class: inline
+    :class: inline matrix
     :header: 0,1,2,3,4,5,6
 
     10.0,2.0,,40.0,-5.0,,24.0
 
 .. csv-table:: w << select.value(v >= 5)
-    :class: inline
+    :class: inline matrix
     :header: 0,1,2,3,4,5,6
 
     10.0,,,40.0,,,24.0
@@ -678,15 +701,16 @@ A monoid or aggregator is used to perform the reduction.
     w << A.reduce_columnwise("times")
 
 .. csv-table:: A
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2,3
+    :stub-columns: 1
 
     **0**,,2.0,,5.0
     **1**,1.5,4.0,,
     **2**,0.5,-7.0,,
 
 .. csv-table:: w << A.reduce_columnwise('times')
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2,3
 
     ,0.75,-56.0,,5.0
@@ -705,15 +729,16 @@ A monoid or aggregator is used to perform the reduction.
     s << A.reduce_scalar("max")
 
 .. csv-table:: A
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2,3
+    :stub-columns: 1
 
     **0**,,2.0,,5.0
     **1**,1.5,4.0,,
     **2**,0.5,-7.0,,
 
 .. csv-table:: s << A.reduce_scalar('max')
-    :class: inline
+    :class: inline matrix
     :header: ,,,,
 
     5.0
@@ -730,13 +755,13 @@ A monoid or aggregator is used to perform the reduction.
     s << gb.agg.argmin(v)
 
 .. csv-table:: v
-    :class: inline
+    :class: inline matrix
     :header: 0,1,2,3,4,5,6
 
     10.0,2.0,,40.0,-5.0,,24.0
 
 .. csv-table:: s << argmin(v)
-    :class: inline
+    :class: inline matrix
     :header: ,,,
 
     4
@@ -761,16 +786,18 @@ To force the transpose to be computed by itself, use it by itself as the right-h
     C << A.T
 
 .. csv-table:: A
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2,3
+    :stub-columns: 1
 
     **0**,,2.0,,5.0
     **1**,1.5,4.0,,
     **2**,0.5,,-7.0,
 
 .. csv-table:: C << A.T
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,,1.5,0.5
     **1**,2.0,4.0,
@@ -802,23 +829,26 @@ The Kronecker product uses a binary operator.
     C << A.kronecker(B, "times")
 
 .. csv-table:: A
-    :class: inline
+    :class: inline matrix
     :header: ,0,1
+    :stub-columns: 1
 
     **0**,1.0,-2.0
     **1**,3.0,
 
 .. csv-table:: B
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2
+    :stub-columns: 1
 
     **0**,,2.0,5.0
     **1**,1.5,4.0,
     **2**,0.5,,-7.0
 
 .. csv-table:: C << A.kronecker(B, 'times')
-    :class: inline
+    :class: inline matrix
     :header: ,0,1,2,3,4,5
+    :stub-columns: 1
 
     **0**,,2.0,5.0,,-4.0,-10.0
     **1**,1.5,4.0,,-3.0,-8.0,
diff --git a/docs/user_guide/operators.rst b/docs/user_guide/operators.rst
index ec28e2fba..8bb5e9fa8 100644
--- a/docs/user_guide/operators.rst
+++ b/docs/user_guide/operators.rst
@@ -314,12 +314,14 @@ each symbol. Each is detailed below.
 The following objects will be used to demonstrate the behavior.
 
 .. csv-table:: Vector v
+    :class: matrix
     :header: 0,1,2,3,4,5
 
     1.0,,2.0,3.5,,9.0
 
 
 .. csv-table:: Vector w
+    :class: matrix
     :header: 0,1,2,3,4,5
 
     7.0,5.2,,3.0,,2.5
@@ -343,6 +345,7 @@ Addition performs an element-wise union between collections, adding overlapping
     v + w
 
 .. csv-table::
+    :class: matrix
     :header: 0,1,2,3,4,5
 
     8.0,5.2,2.0,6.5,,11.5
@@ -358,6 +361,7 @@ and negating any standalone elements from the right-hand object.
     v - w
 
 .. csv-table::
+    :class: matrix
     :header: 0,1,2,3,4,5
 
     -6.0,-5.2,2.0,0.5,,6.5
@@ -373,6 +377,7 @@ overlapping elements.
     v * w
 
 .. csv-table::
+    :class: matrix
     :header: 0,1,2,3,4,5
 
     7.0,,,10.5,,22.5
@@ -392,6 +397,7 @@ elements and always results in a floating-point dtype.
     v / w
 
 .. csv-table::
+    :class: matrix
     :header: 0,1,2,3,4,5
 
     0.142857,,,1.166667,,3.6
@@ -407,6 +413,7 @@ Dividing by zero with floor division will raise a ``ZeroDivisionError``.
     v // w
 
 .. csv-table::
+    :class: matrix
     :header: 0,1,2,3,4,5
 
     0.0,,,1.0,,3.0
@@ -422,6 +429,7 @@ of dividing overlapping elements.
     v % w
 
 .. csv-table::
+    :class: matrix
     :header: 0,1,2,3,4,5
 
     1.0,,,0.5,,1.5
@@ -437,6 +445,7 @@ the power of y for overlapping elements.
     v**w
 
 .. csv-table::
+    :class: matrix
     :header: 0,1,2,3,4,5
 
     1.0,,,42.875,,243.0
@@ -455,6 +464,7 @@ rather than ``all(A == B)``
     v > w
 
 .. csv-table::
+    :class: matrix
     :header: 0,1,2,3,4,5
 
     False,,,True,,True
@@ -464,6 +474,7 @@ rather than ``all(A == B)``
     v == w
 
 .. csv-table::
+    :class: matrix
     :header: 0,1,2,3,4,5
 
     False,,,False,,False
diff --git a/docs/user_guide/udf.rst b/docs/user_guide/udf.rst
index b96097a85..e7b984b44 100644
--- a/docs/user_guide/udf.rst
+++ b/docs/user_guide/udf.rst
@@ -27,6 +27,7 @@ Example user-defined UnaryOp:
     w = v.apply(unary.force_odd).new()
 
 .. csv-table:: w
+    :class: matrix
     :header: 0,1,2,3,4,5
 
     1,3,,3,9,15
@@ -48,6 +49,7 @@ Example lambda usage:
     v.apply(lambda x: x % 5 - 2).new()
 
 .. csv-table::
+    :class: matrix
     :header: 0,1,2,3,4,5
 
     -1,0,,1,1,2