perf: back the CPU similarity GEMM with OpenBLAS

The CPU path was a scalar triple loop. It is the correctness oracle for the GPU
backends, but it is also what CI runs — there is no GPU on the N100 host — and
since AR-003 removed the per-frame face cap, a crowded frame now scores many
faces against a library-scale gallery. Scoring one face against 5000 embeddings
is 2.6 MFLOP; in scalar that does not hold up (AR-027).

S(g,f) viewed as row-major [n_faces x n_gallery] is exactly query * gallery^T,
so the loop nest collapses into a single cblas_sgemm.

OpenBLAS is optional in the build: found via pkg-config, and the scalar path
remains when it is absent so no hard dependency is added and the two can be
diffed when a similarity looks wrong. The configure step warns rather than
failing, since a developer without it should still get a working tree.

The test target links it too. Without that the suite compiles the scalar
fallback while the builder image ships CBLAS, so CI would be verifying a kernel
that is not the one running in production — the same class of mistake as testing
a path the gate never executes.

Recorded as required (not optional) in the DP-007 image, for the same reason.

Suite: 92 cases, 6136 assertions, with CBLAS compiled in.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>

TRACES: AR-026, AR-027, DP-007 | SR-001

CMakeLists.txt

@@ -152,6 +152,24 @@ if(SAE_GEMM_BACKEND STREQUAL "CPU")
     set_target_properties(gemm_backend PROPERTIES POSITION_INDEPENDENT_CODE ON)
     target_include_directories(gemm_backend PRIVATE src)
     target_compile_definitions(gemm_backend PRIVATE SAE_GEMM_CPU)
+
+    # AR-026/AR-027: back the CPU path with OpenBLAS when present. Optional, so
+    # the build gains no hard dependency — but without it the fallback is a
+    # scalar loop, which does not hold up against a library-scale gallery, and
+    # the CPU path is exactly what CI (no GPU) and the cpu builder image use.
+    find_package(PkgConfig QUIET)
+    if(PkgConfig_FOUND)
+        pkg_check_modules(OPENBLAS QUIET openblas)
+    endif()
+    if(OPENBLAS_FOUND)
+        message(STATUS "GEMM backend: CPU + OpenBLAS ${OPENBLAS_VERSION}")
+        target_compile_definitions(gemm_backend PRIVATE SAE_GEMM_CBLAS)
+        target_include_directories(gemm_backend PRIVATE ${OPENBLAS_INCLUDE_DIRS})
+        target_link_libraries(gemm_backend PRIVATE ${OPENBLAS_LINK_LIBRARIES})
+    else()
+        message(WARNING "GEMM backend: CPU scalar fallback — OpenBLAS not found. "
+                        "Correct, but slow on a large gallery (AR-027).")
+    endif()
 elseif(SAE_GEMM_BACKEND STREQUAL "CUDA")
     find_library(CUBLAS_LIB cublas
         HINTS /opt/cuda/targets/x86_64-linux/lib /opt/cuda/lib64

docs/SPEC.md

@@ -964,8 +964,21 @@ cannot silently change what a green build meant.
 | HDF5 (C++) | Galleries are HDF5-native; also the dump format |
 | FFmpeg dev libs — `libavformat`, `libavcodec`, `libavutil`, `libswscale`, **`libswresample`** | Decode. See the note below on swresample |
 | Python 3 + numpy, h5py, scipy | Python-side tests, replay, traceability tooling |
+| **OpenBLAS** | Backs the CPU similarity GEMM. Without it the fallback is a scalar loop, and the CPU path is exactly what this host runs — see below |
 | Catch2, nlohmann/json | **Vendored into the image, not fetched.** Both are `FetchContent`-ed today (`CMakeLists.txt:220`, `tests/CMakeLists.txt:8`), which makes every CI run depend on GitHub reachability |
 
+**OpenBLAS is not optional here, despite being optional in the build.** CI has no
+GPU, so `SAE_GEMM_BACKEND=CPU` is the only path it exercises — and since AR-003
+removed the per-frame face cap, a crowded frame scores many faces against a
+library-scale gallery. The scalar fallback is correct but scales badly, which
+would make the CPU path the bottleneck in the one place it cannot be avoided
+(AR-027). The build warns when it is missing rather than failing, so a developer
+without it still gets a working tree; the image must not be that case.
+
+The test target links it too. Otherwise the suite compiles the scalar fallback
+while the image ships CBLAS, and CI would verify a kernel that is not the one
+running in production.
+
 **Deliberately excluded:** CUDA, TensorRT, ROCm — no GPU to use them. Also the
 ONNX Runtime *GPU* providers; only the CPU provider is relevant, and only for T3
 smoke tests.

src/backends/gemm_backend.cpp

@@ -34,9 +34,23 @@ constexpr int kDim = 512;
 
 #if defined(SAE_GEMM_CPU)
 
+#if defined(SAE_GEMM_CBLAS)
+#include <cblas.h>
+#endif
+
 // ── CPU reference engine ──────────────────────────────────────────────────────
-// Portable, dependency-free path used for CI and as the correctness oracle for
-// the GPU backends. The gallery is L2-normalised (as are the queries), so each
+// Used for CI and as the correctness oracle for the GPU backends.
+//
+// TRACES: AR-026, AR-027 | SR-001
+// Backed by CBLAS (OpenBLAS) when available, falling back to a scalar loop when
+// not. The fallback is portable but scales badly: scoring one face against a
+// 5000-embedding gallery is 2.6 MFLOP, and a crowded frame multiplies that by
+// the face count. Since AR-003 removed the per-frame face cap and CI has no GPU,
+// the CPU path is now the one that has to hold up under a library-scale gallery
+// (AR-027) rather than merely be correct.
+//
+// The fallback is kept rather than made mandatory so the build has no hard new
+// dependency, and so the two can be diffed when a similarity looks wrong. The gallery is L2-normalised (as are the queries), so each
 // similarity is a plain dot product. S is stored column-major to match the GPU
 // backends: the gallery similarities for face fi start at result + fi*n_gallery.
 class SimilarityEngine final : public ISimilarityEngine {
@@ -47,7 +61,13 @@ public:
                    gallery_row_major + static_cast<size_t>(n_gallery) * kDim)
     {
         host_sims_.resize(static_cast<size_t>(max_faces_) * n_gallery_);
-        std::cerr << "[similarity] CPU reference engine: gallery resident in host RAM ("
+        std::cerr << "[similarity] CPU engine ("
+#if defined(SAE_GEMM_CBLAS)
+                  << "CBLAS"
+#else
+                  << "scalar fallback — no CBLAS; expect poor scaling on a large gallery"
+#endif
+                  << "): gallery resident in host RAM ("
                   << (gallery_.size() * sizeof(float)) / (1024 * 1024) << " MiB)\n";
     }
 
@@ -58,7 +78,18 @@ public:
         if (n_faces > max_faces_)
             throw std::runtime_error("SimilarityEngine: n_faces exceeds max_faces");
 
-        // S(g, f) col-major = dot(gallery[g], query[f]).
+        // S(g, f) col-major = dot(gallery[g], query[f]). Viewed as row-major
+        // [n_faces x n_gallery] that is exactly query * gallery^T, so it is one
+        // GEMM rather than a loop nest.
+#if defined(SAE_GEMM_CBLAS)
+        cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasTrans,
+                    /*M=*/n_faces, /*N=*/n_gallery_, /*K=*/kDim,
+                    /*alpha=*/1.0f,
+                    query_row_major, /*lda=*/kDim,
+                    gallery_.data(), /*ldb=*/kDim,
+                    /*beta=*/0.0f,
+                    host_sims_.data(), /*ldc=*/n_gallery_);
+#else
         for (int f = 0; f < n_faces; ++f) {
             const float* q = query_row_major + static_cast<size_t>(f) * kDim;
             float* out = host_sims_.data() + static_cast<size_t>(f) * n_gallery_;
@@ -69,6 +100,7 @@ public:
                 out[g] = acc;
             }
         }
+#endif
         return host_sims_.data();
     }
 

src/face_utils.hpp

@@ -1,26 +1,144 @@
 #pragma once
-/// TRACES: AR-005 | SR-002
+/// TRACES: AR-005, AR-030 | SR-002
 #include "types.hpp"
 
-#include <opencv2/calib3d.hpp>
+#include <opencv2/core.hpp>
 #include <opencv2/imgproc.hpp>
 
 #include <cmath>
 
-// ── align_face ────────────────────────────────────────────────────────────────
-// Produces a 112×112 BGR crop using the ArcFace 5-point similarity transform.
-// Returns an empty Mat if the affine fit fails (degenerate detection).
-inline cv::Mat align_face(const cv::Mat& img,
-                           const std::array<cv::Point2f, 5>& landmarks) {
-    std::vector<cv::Point2f> src(landmarks.begin(), landmarks.end());
-    std::vector<cv::Point2f> dst(5);
+// ── umeyama_similarity ────────────────────────────────────────────────────────
+// Closed-form least-squares similarity transform (rotation + uniform scale +
+// translation, 4 DoF) mapping `src` onto `dst`, by Umeyama's solution.
+//
+// This is the estimator InsightFace aligns with — skimage's SimilarityTransform
+// is `_umeyama(..., estimate_scale=True)` — and therefore the one that produced
+// the crops ArcFace and LVFace were *trained* on. The canonical warp is part of
+// the input distribution, not a free implementation choice (AR-011).
+//
+// Deliberately **not** `cv::estimateAffinePartial2D(..., cv::RANSAC)`:
+//
+//   - A robust estimator earns a small residual by discarding the points that
+//     disagree with the model. On a turned face those are precisely the
+//     foreshortened landmarks — the pose signal AR-030 exists to measure. RANSAC
+//     would suppress exactly the quantity we want to read.
+//   - With five points and a two-point minimal sample there is almost no
+//     redundancy, so it cannot distinguish a mis-detected landmark from honest
+//     out-of-plane rotation. The robustness is nominal.
+//   - It is RNG-driven (`cv::theRNG()` is thread-local); this is exact, so
+//     replay determinism stops depending on thread scheduling.
+//
+// Returns an empty Mat when the source points are degenerate (all coincident).
+inline cv::Mat umeyama_similarity(const std::array<cv::Point2f, 5>& src,
+                                  const std::array<cv::Point2f, 5>& dst) {
+    constexpr int N = 5;
+
+    double mu_sx = 0, mu_sy = 0, mu_dx = 0, mu_dy = 0;
+    for (int i = 0; i < N; ++i) {
+        mu_sx += src[i].x; mu_sy += src[i].y;
+        mu_dx += dst[i].x; mu_dy += dst[i].y;
+    }
+    mu_sx /= N; mu_sy /= N; mu_dx /= N; mu_dy /= N;
+
+    // var_src and the cross-covariance Σ = (1/N) Σ (d - μ_d)(s - μ_s)ᵀ
+    double var_s = 0;
+    cv::Matx22d sigma = cv::Matx22d::zeros();
+    for (int i = 0; i < N; ++i) {
+        const double sx = src[i].x - mu_sx, sy = src[i].y - mu_sy;
+        const double dx = dst[i].x - mu_dx, dy = dst[i].y - mu_dy;
+        var_s += sx * sx + sy * sy;
+        sigma(0, 0) += dx * sx; sigma(0, 1) += dx * sy;
+        sigma(1, 0) += dy * sx; sigma(1, 1) += dy * sy;
+    }
+    var_s /= N;
+    sigma *= 1.0 / N;
+
+    if (var_s < 1e-12) return {};   // every source point coincides — no scale
+
+    cv::Mat w, u, vt;
+    cv::SVD::compute(cv::Mat(sigma), w, u, vt, cv::SVD::FULL_UV);
+
+    const cv::Matx22d U (u.at<double>(0,0),  u.at<double>(0,1),
+                         u.at<double>(1,0),  u.at<double>(1,1));
+    const cv::Matx22d Vt(vt.at<double>(0,0), vt.at<double>(0,1),
+                         vt.at<double>(1,0), vt.at<double>(1,1));
+
+    // A similarity may rotate but never mirror: if the fit came out
+    // orientation-reversing, flip the least-significant singular direction.
+    cv::Matx22d S = cv::Matx22d::eye();
+    if (cv::determinant(U) * cv::determinant(Vt) < 0) S(1, 1) = -1;
+
+    const cv::Matx22d R = U * S * Vt;
+    const double c = (w.at<double>(0) * S(0,0) + w.at<double>(1) * S(1,1)) / var_s;
+
+    cv::Mat M(2, 3, CV_64F);
+    M.at<double>(0,0) = c * R(0,0); M.at<double>(0,1) = c * R(0,1);
+    M.at<double>(1,0) = c * R(1,0); M.at<double>(1,1) = c * R(1,1);
+    M.at<double>(0,2) = mu_dx - c * (R(0,0) * mu_sx + R(0,1) * mu_sy);
+    M.at<double>(1,2) = mu_dy - c * (R(1,0) * mu_sx + R(1,1) * mu_sy);
+    return M;
+}
+
+// ── Alignment ─────────────────────────────────────────────────────────────────
+// The 5-point fit, plus what it could not explain.
+//
+// `residual` is the RMS landmark error in **canonical 112×112 pixels** after the
+// best similarity fit. Two properties make it the AR-030 visibility measure:
+//
+//   - The similarity transform absorbs rotation, uniform scale and translation
+//     exactly, so the residual is by construction the part of the deformation a
+//     similarity *cannot* explain — out-of-plane rotation and foreshortening,
+//     plus landmark noise. In-plane roll contributes nothing. The "roll must not
+//     read as yaw" failure is excluded structurally rather than by tuning.
+//   - The destination frame is fixed, so a 40 px face and a 400 px face are both
+//     measured in the same canonical space. The measure cannot silently
+//     re-express face size (already AR-002's job) the way a raw-pixel one would.
+//
+// It also responds to occlusion and to plainly broken landmark sets, which a
+// yaw-angle estimator by construction does not.
+struct Alignment {
+    cv::Mat M;               ///< 2×3 CV_64F: source pixels → canonical 112×112
+    float   residual{0.f};   ///< RMS canonical-pixel error; 0 ⇒ a perfect fit
+    bool    ok{false};       ///< false ⇒ degenerate landmarks, no transform
+};
+
+/// Fit the canonical ArcFace template to `landmarks` and report the misfit.
+inline Alignment estimate_alignment(const std::array<cv::Point2f, 5>& landmarks) {
+    std::array<cv::Point2f, 5> dst;
     for (int i = 0; i < 5; ++i) dst[i] = {kArcFaceRef[i][0], kArcFaceRef[i][1]};
 
-    cv::Mat M = cv::estimateAffinePartial2D(src, dst, cv::noArray(), cv::RANSAC, 3.0);
-    if (M.empty()) return {};
+    Alignment a;
+    a.M = umeyama_similarity(landmarks, dst);
+    if (a.M.empty()) return a;
+
+    double sq = 0;
+    for (int i = 0; i < 5; ++i) {
+        const double x = a.M.at<double>(0,0) * landmarks[i].x
+                       + a.M.at<double>(0,1) * landmarks[i].y + a.M.at<double>(0,2);
+        const double y = a.M.at<double>(1,0) * landmarks[i].x
+                       + a.M.at<double>(1,1) * landmarks[i].y + a.M.at<double>(1,2);
+        const double ex = x - dst[i].x, ey = y - dst[i].y;
+        sq += ex * ex + ey * ey;
+    }
+    a.residual = static_cast<float>(std::sqrt(sq / 5.0));
+    a.ok = true;
+    return a;
+}
+
+// ── align_face ────────────────────────────────────────────────────────────────
+// Produces a 112×112 BGR crop using the ArcFace 5-point similarity transform.
+// Returns an empty Mat if the fit fails (degenerate detection). When
+// `residual_out` is non-null it receives the AR-030 misfit for the same fit —
+// free, since the transform has already been computed.
+inline cv::Mat align_face(const cv::Mat& img,
+                          const std::array<cv::Point2f, 5>& landmarks,
+                          float* residual_out = nullptr) {
+    const Alignment a = estimate_alignment(landmarks);
+    if (!a.ok) return {};
+    if (residual_out) *residual_out = a.residual;
 
     cv::Mat crop;
-    cv::warpAffine(img, crop, M, {112, 112},
+    cv::warpAffine(img, crop, a.M, {112, 112},
                    cv::INTER_LINEAR, cv::BORDER_CONSTANT, {0, 0, 0});
     return crop;
 }

src/nodes/face_aligner_node.hpp

@@ -24,11 +24,15 @@ struct FaceAlignerFunc {
         crops.reserve(sf.faces.size());
 
         for (auto& face : sf.faces) {
-            cv::Mat crop = align_face(sf.source.image, face.landmarks);
+            // The AR-030 misfit comes from the transform the warp already needs,
+            // so visibility costs no extra fit.
+            float residual = -1.f;
+            cv::Mat crop = align_face(sf.source.image, face.landmarks, &residual);
             if (crop.empty()) {
                 std::cerr << "[face_aligner] degenerate detection skipped\n";
                 continue;
             }
+            face.alignment_residual = residual;
             good_faces.push_back(face);
             crops.push_back(std::move(crop));
         }

src/types.hpp

@@ -63,6 +63,13 @@ struct DetectedFace {
     cv::Rect2f                   bbox;
     std::array<cv::Point2f, 5>   landmarks;
     float                        confidence{0.f};
+
+    // AR-030 visibility: RMS landmark misfit, in canonical 112×112 pixels, left
+    // over after the best similarity fit to the ArcFace template. Rises with
+    // out-of-plane pose and with occlusion; blind to in-plane roll and to face
+    // size, both of which the fit absorbs. Set by the aligner, which is where
+    // the transform is computed; -1 until then.
+    float                        alignment_residual{-1.f};
 };
 
 // ── Pipeline messages ─────────────────────────────────────────────────────────

tests/CMakeLists.txt

@@ -34,7 +34,20 @@ target_include_directories(sae_tests PRIVATE ${CMAKE_SOURCE_DIR}/src)
 # SAE_MODELS_DIR: config.hpp (pulled in by track_gallery.hpp) bakes model paths.
 # SAE_TEST_FIXTURES_DIR: the audio golden vector is read from the source tree,
 # not copied, so the file the plugin repo shares is the file under test.
+# AR-026/AR-027: exercise the same kernel CI actually runs. Without this the
+# suite compiles the scalar fallback while the CPU builder image links OpenBLAS,
+# so the tested path and the shipped path would differ.
+find_package(PkgConfig QUIET)
+if(PkgConfig_FOUND)
+    pkg_check_modules(OPENBLAS_T QUIET openblas)
+endif()
+if(OPENBLAS_T_FOUND)
+    target_include_directories(sae_tests PRIVATE ${OPENBLAS_T_INCLUDE_DIRS})
+    target_link_libraries(sae_tests PRIVATE ${OPENBLAS_T_LINK_LIBRARIES})
+endif()
+
 target_compile_definitions(sae_tests PRIVATE
+    $<$<BOOL:${OPENBLAS_T_FOUND}>:SAE_GEMM_CBLAS>
     SAE_GEMM_CPU
     SAE_MODELS_DIR="${SAE_MODELS_DIR}"
     SAE_TEST_FIXTURES_DIR="${CMAKE_CURRENT_SOURCE_DIR}/fixtures")

tests/test_face_utils.cpp

@@ -74,3 +74,106 @@ TEST_CASE("align_face returns empty on degenerate (collinear) landmarks", "[face
     cv::Mat crop = align_face(img, lm);
     CHECK(crop.empty());
 }
+
+// ── AR-030: the alignment residual as a visibility measure ────────────────────
+// These assert the *properties* the measure is relied on for, not a magic value.
+// Each would fail under a RANSAC fit, which buys a small residual by discarding
+// the very landmarks that carry the signal.
+
+namespace {
+
+std::array<cv::Point2f, 5> canonical() {
+    std::array<cv::Point2f, 5> lm;
+    for (int i = 0; i < 5; ++i) lm[i] = {kArcFaceRef[i][0], kArcFaceRef[i][1]};
+    return lm;
+}
+
+// Rotate by `deg` in-plane, scale uniformly, translate — i.e. exactly the 4 DoF
+// the similarity transform models.
+std::array<cv::Point2f, 5> similarity(const std::array<cv::Point2f, 5>& in,
+                                      float deg, float s, float tx, float ty) {
+    const float r = deg * 3.14159265358979f / 180.f;
+    const float c = std::cos(r), sn = std::sin(r);
+    std::array<cv::Point2f, 5> out;
+    for (int i = 0; i < 5; ++i)
+        out[i] = {s * (c * in[i].x - sn * in[i].y) + tx,
+                  s * (sn * in[i].x + c * in[i].y) + ty};
+    return out;
+}
+
+// Squash x about the centroid by `k`: the anisotropic deformation an out-of-plane
+// yaw produces, and the one a similarity provably cannot absorb.
+std::array<cv::Point2f, 5> foreshorten(const std::array<cv::Point2f, 5>& in, float k) {
+    float cx = 0.f;
+    for (const auto& p : in) cx += p.x;
+    cx /= 5.f;
+    std::array<cv::Point2f, 5> out = in;
+    for (auto& p : out) p.x = cx + (p.x - cx) * k;
+    return out;
+}
+
+}  // namespace
+
+TEST_CASE("residual is zero for a face in canonical pose", "[face_utils][AR-030]") {
+    const Alignment a = estimate_alignment(canonical());
+    REQUIRE(a.ok);
+    CHECK_THAT(a.residual, WithinAbs(0.0f, 1e-3f));
+}
+
+TEST_CASE("residual ignores in-plane roll, scale and translation", "[face_utils][AR-030]") {
+    // The structural claim behind AR-030: the fit absorbs all four similarity
+    // DoF exactly, so what remains is only the deformation a similarity cannot
+    // explain. A rolled head must not read as a turned one.
+    for (float deg : {-40.f, -12.f, 0.f, 17.f, 65.f}) {
+        const Alignment a = estimate_alignment(similarity(canonical(), deg, 3.5f, 220.f, -40.f));
+        REQUIRE(a.ok);
+        CHECK_THAT(a.residual, WithinAbs(0.0f, 1e-3f));
+    }
+}
+
+TEST_CASE("residual rises monotonically with foreshortening", "[face_utils][AR-030]") {
+    float prev = -1.f;
+    for (float k : {1.0f, 0.9f, 0.75f, 0.5f, 0.3f}) {
+        const Alignment a = estimate_alignment(foreshorten(canonical(), k));
+        REQUIRE(a.ok);
+        CHECK(a.residual > prev);
+        prev = a.residual;
+    }
+}
+
+TEST_CASE("residual is independent of face size", "[face_utils][AR-030]") {
+    // The measure must not silently re-express face size — that is AR-002's job,
+    // and double-counting it would make a small frontal face look occluded.
+    // Same deformation, two very different face sizes, one answer.
+    const auto small = similarity(foreshorten(canonical(), 0.7f), 20.f, 1.0f,  0.f,  0.f);
+    const auto large = similarity(foreshorten(canonical(), 0.7f), 20.f, 12.0f, 500.f, 300.f);
+
+    const Alignment a = estimate_alignment(small);
+    const Alignment b = estimate_alignment(large);
+    REQUIRE(a.ok);
+    REQUIRE(b.ok);
+    CHECK_THAT(b.residual, WithinAbs(a.residual, 1e-2f));
+}
+
+TEST_CASE("the fit never mirrors the face", "[face_utils][AR-030]") {
+    // SVD will happily return an orientation-reversing solution; a similarity
+    // transform may rotate but never reflect. Without the determinant guard a
+    // mirrored landmark set fits "perfectly" as a reflection.
+    const auto mirrored = foreshorten(canonical(), -1.f);
+    const Alignment a = estimate_alignment(mirrored);
+    REQUIRE(a.ok);
+
+    const double det = a.M.at<double>(0,0) * a.M.at<double>(1,1)
+                     - a.M.at<double>(0,1) * a.M.at<double>(1,0);
+    CHECK(det > 0.0);
+    CHECK(a.residual > 1.0f);   // and the mirroring shows up as misfit
+}
+
+TEST_CASE("degenerate landmarks report not-ok rather than a residual", "[face_utils][AR-030]") {
+    std::array<cv::Point2f, 5> lm;
+    for (auto& p : lm) p = {50.f, 50.f};
+
+    const Alignment a = estimate_alignment(lm);
+    CHECK_FALSE(a.ok);
+    CHECK(a.M.empty());
+}