diff --git a/benchmarks/broadcast_small_inv.cpp b/benchmarks/broadcast_small_inv.cpp
new file mode 100644
index 000000000..295ca1426
--- /dev/null
+++ b/benchmarks/broadcast_small_inv.cpp
@@ -0,0 +1,24 @@
+#include "./bench_common.hpp"
+
+template <int N, long dim>
+static void inv(benchmark::State &state) {
+  nda::array<double, 3> W(dim, N, N), Wi(dim, N, N);
+  for (int k = 0; k < dim; ++k) {
+      for (int i = 0; i < N; ++i)
+        for (int j = 0; j < N; ++j) W(k, i, j) = (i > j ? 0.5 + i + 2.5 * j : i * 0.8 - j - 0.5);
+  }
+
+  while (state.KeepRunning()) {
+    Wi = inverse(W);
+  }
+}
+
+BENCHMARK_TEMPLATE(inv, 1, 100);
+BENCHMARK_TEMPLATE(inv, 1, 10000);
+BENCHMARK_TEMPLATE(inv, 1, 1000000);
+BENCHMARK_TEMPLATE(inv, 2, 100);
+BENCHMARK_TEMPLATE(inv, 2, 10000);
+BENCHMARK_TEMPLATE(inv, 2, 1000000);
+BENCHMARK_TEMPLATE(inv, 3, 100);
+BENCHMARK_TEMPLATE(inv, 3, 10000);
+BENCHMARK_TEMPLATE(inv, 3, 1000000);
diff --git a/c++/nda/linalg/det_and_inverse.hpp b/c++/nda/linalg/det_and_inverse.hpp
index a804351fa..659067c0c 100644
--- a/c++/nda/linalg/det_and_inverse.hpp
+++ b/c++/nda/linalg/det_and_inverse.hpp
@@ -173,6 +173,14 @@ namespace nda {
     return r;
   }
 
+  template <Array A>
+  auto inverse(A const &a) {
+      auto r = make_regular(a);
+      auto long_axis = stdutil::mpop<2>(r.indexmap().lengths());
+      for_each(long_axis, [&r, _ = range()](auto &&... i) {inverse_in_place(make_matrix_view(r(i..., _, _))); });
+      return r;
+  }
+
 } // namespace nda
 
 namespace nda::clef {
diff --git a/test/c++/nda_linear_algebra.cpp b/test/c++/nda_linear_algebra.cpp
index acdbea2ee..ea4174ebe 100644
--- a/test/c++/nda_linear_algebra.cpp
+++ b/test/c++/nda_linear_algebra.cpp
@@ -327,6 +327,30 @@ TEST(Inverse, Small) { //NOLINT
 
 // ==============================================================
 
+TEST(Inverse, BroadcastInverse) { //NOLINT
+
+  for (auto N : {1, 2, 3, 4}) {
+
+    long dim = 100000;
+    nda::array<double, 3> W(dim, N, N);
+
+    for (int k = 0; k < dim; ++k) {
+        for (int i = 0; i < N; ++i)
+            for (int j = 0; j < N; ++j) W(k, i, j) = (i > j ? 0.5 + i + 2.5 * j : i * 0.8 - j - 0.5);
+      }
+
+    auto Wi = inverse(W);
+    auto Wii = inverse(Wi);
+    for ( int k = 0; k < dim; ++k) {
+        EXPECT_NEAR(determinant(make_matrix_view(Wi(k, range(), range()))), 1.0/determinant(make_matrix_view(W(k, range(), range()))), 1.e-12);
+        EXPECT_ARRAY_NEAR(make_matrix_view(W(k,range(),range())) * make_matrix_view(Wi(k, range(), range())), nda::eye<double>(N), 1.e-13);
+        EXPECT_ARRAY_NEAR(make_matrix_view(Wii(k,range(),range())), make_matrix_view(W(k, range(), range())), 1.e-12);
+    }
+  }
+}
+
+// ==============================================================
+
 TEST(Matvecmul, Promotion) { //NOLINT
 
   matrix<int> Ai   = {{1, 2}, {3, 4}};