benchmark the `linalg` module (dragoon/komodo!43)

this MR
- adds `criterion` as a dependency
- creates a `linalg.rs` benchmark file
- makes the following function `pub`lic
  - `Matrix::transpose`
  - `Matrix::invert`
  - `Matrix::mul`
- creates a new `benches/` directory containing
  - a README with commands
  - a `plot.py` file to plot results
  - a `linalg.rs` file with the benchmarks

## example results
![Figure_1](/uploads/f352a6f411662361fa9ca381710271d5/Figure_1.png)

.gitignore

 target/
 Cargo.lock
+
+*.ndjson

Cargo.toml

@@ -26,3 +26,7 @@ criterion = "0.3"
 [[bench]]
 name = "recoding"
 harness = false
+
+[[bench]]
+name = "linalg"
+harness = false

benches/README.md

+## run the benchmarks
+```shell
+nushell> cargo criterion --output-format verbose --message-format json out> results.ndjson
+```
+
+## plot the results
+```shell
+python benches/plot.py results.ndjson
+```

benches/linalg.rs

+use ark_bls12_381::Bls12_381;
+use ark_ec::pairing::Pairing;
+use criterion::{black_box, criterion_group, criterion_main, Criterion};
+
+use komodo::linalg::Matrix;
+
+fn inverse_template<E: Pairing>(c: &mut Criterion, n: usize) {
+    let matrix = Matrix::<E::ScalarField>::random(n, n);
+
+    c.bench_function(
+        &format!("inverse {}x{} on {}", n, n, std::any::type_name::<E>()),
+        |b| b.iter(|| matrix.invert().unwrap()),
+    );
+}
+
+fn inverse(c: &mut Criterion) {
+    for n in [10, 15, 20, 30, 40, 60, 80, 120, 160, 240, 320] {
+        inverse_template::<Bls12_381>(c, black_box(n));
+    }
+}
+
+fn transpose_template<E: Pairing>(c: &mut Criterion, n: usize) {
+    let matrix = Matrix::<E::ScalarField>::random(n, n);
+
+    c.bench_function(
+        &format!("transpose {}x{} on {}", n, n, std::any::type_name::<E>()),
+        |b| b.iter(|| matrix.transpose()),
+    );
+}
+
+fn transpose(c: &mut Criterion) {
+    for n in [10, 15, 20, 30, 40, 60, 80, 120, 160, 240, 320] {
+        transpose_template::<Bls12_381>(c, black_box(n));
+    }
+}
+
+fn mul_template<E: Pairing>(c: &mut Criterion, n: usize) {
+    let mat_a = Matrix::<E::ScalarField>::random(n, n);
+    let mat_b = Matrix::<E::ScalarField>::random(n, n);
+
+    c.bench_function(
+        &format!("mul {}x{} on {}", n, n, std::any::type_name::<E>()),
+        |b| b.iter(|| mat_a.mul(&mat_b)),
+    );
+}
+
+fn mul(c: &mut Criterion) {
+    for n in [10, 15, 20, 30, 40, 60, 80, 120, 160, 240, 320] {
+        mul_template::<Bls12_381>(c, black_box(n));
+    }
+}
+
+criterion_group!(benches, inverse, transpose, mul);
+criterion_main!(benches);

benches/plot.py

+import matplotlib.pyplot as plt
+import json
+import sys
+import os
+from typing import Any, Dict, List
+
+NB_NS_IN_MS = 1e6
+
+Data = List[Dict[str, Any]]
+
+
+def extract(data: Data, k1: str, k2: str) -> List[float]:
+    return [line[k1][k2] / NB_NS_IN_MS for line in data]
+
+
+def plot(data: Data, key: str, ax):
+    filtered_data = list(filter(lambda line: line["id"].startswith(key), data))
+
+    sizes = [
+        int(line["id"].split(' ')[1].split('x')[0]) for line in filtered_data
+    ]
+
+    means = extract(filtered_data, "mean", "estimate")
+    up = extract(filtered_data, "mean", "upper_bound")
+    down = extract(filtered_data, "mean", "lower_bound")
+
+    ax.plot(sizes, means, label="mean", color="blue")
+    ax.fill_between(sizes, down, up, color="blue", alpha=0.3, label="mean bounds")
+
+    medians = extract(filtered_data, "median", "estimate")
+    up = extract(filtered_data, "median", "upper_bound")
+    down = extract(filtered_data, "median", "lower_bound")
+
+    ax.plot(sizes, medians, label="median", color="orange")
+    ax.fill_between(sizes, down, up, color="orange", alpha=0.3, label="median bounds")
+
+
+def parse_args():
+    if len(sys.argv) == 1:
+        print("please give a filename as first positional argument")
+        exit(1)
+
+    return sys.argv[1]
+
+
+def read_data(data_file: str) -> Data:
+    if not os.path.exists(data_file):
+        print(f"no such file: `{data_file}`")
+        exit(1)
+
+    with open(data_file, "r") as file:
+        data = list(filter(
+            lambda line: line["reason"] == "benchmark-complete",
+            map(
+                json.loads,
+                file.readlines()
+            )
+        ))
+
+    return data
+
+
+if __name__ == "__main__":
+    results_file = parse_args()
+    data = read_data(results_file)
+
+    labels = ["transpose", "mul", "inverse"]
+
+    fig, axs = plt.subplots(len(labels), 1)
+
+    for label, ax in zip(labels, axs):
+        plot(data, key=label, ax=ax)
+        ax.set_title(label)
+        ax.set_yscale("log")
+        ax.set_ylabel("time in ms")
+        ax.legend()
+        ax.grid()
+
+    plt.show()

src/linalg.rs

@@ -135,7 +135,7 @@ impl<T: Field> Matrix<T> {
         }
     }
 
-    pub(super) fn invert(&self) -> Result<Self, KomodoError> {
+    pub fn invert(&self) -> Result<Self, KomodoError> {
         if self.height != self.width {
             return Err(KomodoError::NonSquareMatrix(self.height, self.width));
         }
@@ -209,7 +209,7 @@ impl<T: Field> Matrix<T> {
         nb_non_zero_rows
     }
 
-    pub(super) fn mul(&self, rhs: &Self) -> Result<Self, KomodoError> {
+    pub fn mul(&self, rhs: &Self) -> Result<Self, KomodoError> {
         if self.width != rhs.height {
             return Err(KomodoError::IncompatibleMatrixShapes(
                 self.height,
@@ -239,7 +239,7 @@ impl<T: Field> Matrix<T> {
         })
     }
 
-    pub(super) fn transpose(&self) -> Self {
+    pub fn transpose(&self) -> Self {
         let height = self.width;
         let width = self.height;