diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h
index 5271643..f977a29 100644
--- a/Eigen/src/Householder/Householder.h
+++ b/Eigen/src/Householder/Householder.h
@@ -236,6 +236,47 @@
   tau = conj((beta - c0) / beta);
 }
 
+namespace internal {
+
+template <typename Derived, typename EssentialPart,
+          bool Fused = !Derived::IsRowMajor && EssentialPart::ColsAtCompileTime == 1 &&
+                       (EssentialPart::RowsAtCompileTime == 1 || EssentialPart::RowsAtCompileTime == 2)>
+struct householder_apply_left_impl {
+  using Scalar = typename Derived::Scalar;
+  static EIGEN_DEVICE_FUNC void run(MatrixBase<Derived>& mat, const EssentialPart& essential, const Scalar& tau,
+                                    Scalar* workspace) {
+    Map<typename plain_row_type<typename Derived::PlainObject>::type> tmp(workspace, mat.cols());
+    Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(mat.derived(), 1, 0,
+                                                                                        mat.rows() - 1, mat.cols());
+    tmp.noalias() = essential.adjoint() * bottom.unwind();
+    tmp = tau * (tmp + mat.row(0));
+    mat.row(0) -= tmp;
+    bottom.unwind().noalias() -= essential * tmp;
+  }
+};
+
+// Two- and three-element reflectors on column-major storage: finishing each column before advancing replaces
+// four strided row passes, which thrash the cache at large outer strides (issue #3160).
+template <typename Derived, typename EssentialPart>
+struct householder_apply_left_impl<Derived, EssentialPart, true> {
+  using Scalar = typename Derived::Scalar;
+  static EIGEN_DEVICE_FUNC void run(MatrixBase<Derived>& mat, const EssentialPart& essential, const Scalar& tau,
+                                    Scalar*) {
+    // Evaluated once so that the column loop reads plain coefficients; tau may reference a coefficient of mat.
+    const Matrix<Scalar, EssentialPart::RowsAtCompileTime, 1> v = essential;
+    const Scalar tauValue = tau;
+    Block<Derived, EssentialPart::RowsAtCompileTime, Derived::ColsAtCompileTime> bottom(mat.derived(), 1, 0,
+                                                                                        mat.rows() - 1, mat.cols());
+    for (Index j = 0; j < mat.cols(); ++j) {
+      const Scalar tmp = tauValue * (v.dot(bottom.col(j)) + mat.coeff(0, j));
+      mat.coeffRef(0, j) -= tmp;
+      bottom.col(j) -= v * tmp;
+    }
+  }
+};
+
+}  // namespace internal
+
 /** Apply the elementary reflector H given by
  * \f$ H = I - tau v v^*\f$
  * with
@@ -258,13 +299,7 @@
   if (rows() == 1) {
     *this *= Scalar(1) - tau;
   } else if (!numext::is_exactly_zero(tau)) {
-    Map<typename internal::plain_row_type<PlainObject>::type> tmp(workspace, cols());
-    Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(derived(), 1, 0, rows() - 1,
-                                                                                        cols());
-    tmp.noalias() = essential.adjoint() * bottom.unwind();
-    tmp = tau * (tmp + this->row(0));
-    this->row(0) = this->row(0) - tmp;
-    bottom.unwind().noalias() -= essential * tmp;
+    internal::householder_apply_left_impl<Derived, EssentialPart>::run(*this, essential, tau, workspace);
   }
 }
 
diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h
index ff92663..06b5978 100644
--- a/Eigen/src/LU/PartialPivLU.h
+++ b/Eigen/src/LU/PartialPivLU.h
@@ -363,6 +363,17 @@
   using BlockType = Ref<Matrix<Scalar, Dynamic, Dynamic, StorageOrder>>;
   using RealScalar = typename MatrixType::RealScalar;
 
+  static void apply_row_transpositions(BlockType& matrix, Index first, Index count, const PivIndex* transpositions) {
+    EIGEN_IF_CONSTEXPR (StorageOrder == ColMajor) {
+      // Keep the pivot rows of one column in cache, even when the outer stride maps every column to the same set.
+      for (Index j = 0; j < matrix.cols(); ++j)
+        for (Index i = first; i < first + count; ++i)
+          numext::swap(matrix.coeffRef(i, j), matrix.coeffRef(transpositions[i], j));
+    } else {
+      for (Index i = first; i < first + count; ++i) matrix.row(i).swap(matrix.row(transpositions[i]));
+    }
+  }
+
   /** \internal performs the LU decomposition in-place of the matrix \a lu
    * using an unblocked algorithm.
    *
@@ -489,15 +500,13 @@
       // update permutations and apply them to A_0
       if (k > 0) {
         BlockType A_0 = lu.block(0, 0, rows, k);
-        for (Index i = k; i < k + bs; ++i) {
-          Index piv = (row_transpositions[i] += internal::convert_index<PivIndex>(k));
-          A_0.row(i).swap(A_0.row(piv));
-        }
+        for (Index i = k; i < k + bs; ++i) row_transpositions[i] += internal::convert_index<PivIndex>(k);
+        apply_row_transpositions(A_0, k, bs, row_transpositions);
       }
 
       if (trows) {
         // apply permutations to A_2
-        for (Index i = k; i < k + bs; ++i) A_2.row(i).swap(A_2.row(row_transpositions[i]));
+        apply_row_transpositions(A_2, k, bs, row_transpositions);
 
         // A12 = A11^-1 A12
         A11.template triangularView<UnitLower>().solveInPlace(A12);
diff --git a/benchmarks/Eigenvalues/CMakeLists.txt b/benchmarks/Eigenvalues/CMakeLists.txt
index 986d701..f59b368 100644
--- a/benchmarks/Eigenvalues/CMakeLists.txt
+++ b/benchmarks/Eigenvalues/CMakeLists.txt
@@ -3,6 +3,7 @@
 
 eigen_add_benchmark(bench_eigensolver bench_eigensolver.cpp)
 eigen_add_benchmark(bench_eigensolver_double bench_eigensolver.cpp DEFINITIONS SCALAR=double)
+eigen_add_benchmark(bench_eigensolver_stride bench_eigensolver_stride.cpp)
 eigen_add_benchmark(bench_eig33 bench_eig33.cpp)
 eigen_add_benchmark(bench_tridiagonal_bisection bench_tridiagonal_bisection.cpp)
 eigen_add_benchmark(bench_tridiagonal_inverse_iteration bench_tridiagonal_inverse_iteration.cpp)
diff --git a/benchmarks/Eigenvalues/bench_eigensolver_stride.cpp b/benchmarks/Eigenvalues/bench_eigensolver_stride.cpp
new file mode 100644
index 0000000..a73f8df
--- /dev/null
+++ b/benchmarks/Eigenvalues/bench_eigensolver_stride.cpp
@@ -0,0 +1,38 @@
+// SPDX-FileCopyrightText: The Eigen Authors
+// SPDX-License-Identifier: MPL-2.0
+
+#include <benchmark/benchmark.h>
+#include <Eigen/Eigenvalues>
+#include <cstdlib>
+
+using namespace Eigen;
+
+template <typename Scalar>
+static void BM_EigenSolverStride(benchmark::State& state) {
+  using Mat = Matrix<Scalar, Dynamic, Dynamic>;
+  const Index n = state.range(0);
+  const bool computeVectors = state.range(1) != 0;
+  std::srand(1);
+  const Mat a = Mat::Random(n, n);
+  EigenSolver<Mat> solver(a, true);
+  if (solver.info() != Success) {
+    state.SkipWithError("EigenSolver did not converge");
+    return;
+  }
+  const Mat vectors = solver.pseudoEigenvectors();
+  const Mat values = solver.pseudoEigenvalueMatrix();
+  const Scalar residual = (a * vectors - vectors * values).norm();
+  const Scalar bound = Scalar(64 * n) * NumTraits<Scalar>::epsilon() * a.norm() * vectors.norm();
+  if (!(residual <= bound)) {
+    state.SkipWithError("Eigenpair residual failed");
+    return;
+  }
+  for (auto _ : state) {
+    solver.compute(a, computeVectors);
+    benchmark::DoNotOptimize(solver.eigenvalues().data());
+    benchmark::ClobberMemory();
+  }
+}
+
+BENCHMARK_TEMPLATE(BM_EigenSolverStride, float)->ArgsProduct({{32, 128, 500, 512, 768, 1000, 1001, 1024}, {0, 1}});
+BENCHMARK_TEMPLATE(BM_EigenSolverStride, double)->ArgsProduct({{32, 128, 500, 512, 768, 1000, 1001, 1024}, {0, 1}});
diff --git a/benchmarks/Householder/CMakeLists.txt b/benchmarks/Householder/CMakeLists.txt
index 8310c92..bb3f283 100644
--- a/benchmarks/Householder/CMakeLists.txt
+++ b/benchmarks/Householder/CMakeLists.txt
@@ -2,3 +2,4 @@
 # SPDX-License-Identifier: MPL-2.0
 
 eigen_add_benchmark(bench_householder bench_householder.cpp)
+eigen_add_benchmark(bench_householder_short bench_householder_short.cpp)
diff --git a/benchmarks/Householder/bench_householder_short.cpp b/benchmarks/Householder/bench_householder_short.cpp
new file mode 100644
index 0000000..37bef4c
--- /dev/null
+++ b/benchmarks/Householder/bench_householder_short.cpp
@@ -0,0 +1,53 @@
+// SPDX-FileCopyrightText: The Eigen Authors
+// SPDX-License-Identifier: MPL-2.0
+
+#include <benchmark/benchmark.h>
+#include <Eigen/Householder>
+#include <cstdlib>
+
+using namespace Eigen;
+
+template <typename Scalar, int Size, int Side>
+static void BM_ShortHouseholder(benchmark::State& state) {
+  using Mat = Matrix<Scalar, Dynamic, Dynamic>;
+  const Index n = state.range(0);
+  const Index stride = n + state.range(1);
+  std::srand(1);
+  Mat storage = Mat::Random(stride, n);
+  Matrix<Scalar, Dynamic, 1> workspace(n);
+  Matrix<Scalar, Size, 1> vector = Matrix<Scalar, Size, 1>::Random();
+  Matrix<Scalar, Size - 1, 1> essential;
+  Scalar tau, beta;
+  vector.makeHouseholder(essential, tau, beta);
+  vector << Scalar(1), essential;
+  const Matrix<Scalar, Size, Size> h = Matrix<Scalar, Size, Size>::Identity() - tau * vector * vector.adjoint();
+  auto block = storage.block(1, 1, Side == OnTheLeft ? Size : n - 1, Side == OnTheLeft ? n - 1 : Size);
+  const Mat original = block;
+  Mat expected;
+  if (Side == OnTheLeft) {
+    expected = h * original;
+    block.applyHouseholderOnTheLeft(essential, tau, workspace.data());
+  } else {
+    expected = original * h;
+    block.applyHouseholderOnTheRight(essential, tau, workspace.data());
+  }
+  if (!((block - expected).norm() <= Scalar(16 * Size) * NumTraits<Scalar>::epsilon() * original.norm())) {
+    state.SkipWithError("Householder application failed");
+    return;
+  }
+  for (auto _ : state) {
+    if (Side == OnTheLeft)
+      block.applyHouseholderOnTheLeft(essential, tau, workspace.data());
+    else
+      block.applyHouseholderOnTheRight(essential, tau, workspace.data());
+    benchmark::DoNotOptimize(block.data());
+    benchmark::ClobberMemory();
+  }
+}
+
+BENCHMARK_TEMPLATE(BM_ShortHouseholder, float, 2, OnTheLeft)->ArgsProduct({{32, 500, 512, 768, 1000, 1024}, {0, 8}});
+BENCHMARK_TEMPLATE(BM_ShortHouseholder, float, 3, OnTheLeft)->ArgsProduct({{32, 500, 512, 768, 1000, 1024}, {0, 8}});
+BENCHMARK_TEMPLATE(BM_ShortHouseholder, double, 2, OnTheLeft)->ArgsProduct({{32, 500, 512, 768, 1000, 1024}, {0, 8}});
+BENCHMARK_TEMPLATE(BM_ShortHouseholder, double, 3, OnTheLeft)->ArgsProduct({{32, 500, 512, 768, 1000, 1024}, {0, 8}});
+BENCHMARK_TEMPLATE(BM_ShortHouseholder, float, 3, OnTheRight)->ArgsProduct({{32, 500, 512, 768, 1000, 1024}, {0, 8}});
+BENCHMARK_TEMPLATE(BM_ShortHouseholder, double, 3, OnTheRight)->ArgsProduct({{32, 500, 512, 768, 1000, 1024}, {0, 8}});
diff --git a/benchmarks/LU/CMakeLists.txt b/benchmarks/LU/CMakeLists.txt
index b6923e8..45aa594 100644
--- a/benchmarks/LU/CMakeLists.txt
+++ b/benchmarks/LU/CMakeLists.txt
@@ -2,4 +2,5 @@
 # SPDX-License-Identifier: MPL-2.0
 
 eigen_add_benchmark(bench_lu bench_lu.cpp)
+eigen_add_benchmark(bench_lu_stride bench_lu_stride.cpp)
 eigen_add_benchmark(bench_rcond bench_rcond.cpp)
diff --git a/benchmarks/LU/bench_lu_stride.cpp b/benchmarks/LU/bench_lu_stride.cpp
new file mode 100644
index 0000000..b72ef59
--- /dev/null
+++ b/benchmarks/LU/bench_lu_stride.cpp
@@ -0,0 +1,44 @@
+// SPDX-FileCopyrightText: The Eigen Authors
+// SPDX-License-Identifier: MPL-2.0
+
+#include <benchmark/benchmark.h>
+#include <Eigen/LU>
+#include <cstdlib>
+
+using namespace Eigen;
+
+template <typename Scalar, int StorageOrder>
+static void BM_PartialPivLUStride(benchmark::State& state) {
+  using Mat = Matrix<Scalar, Dynamic, Dynamic, StorageOrder>;
+  const Index n = state.range(0);
+  const Index stride = n + state.range(1);
+  std::srand(1);
+  const Mat a = Mat::Random(n, n);
+  Matrix<Scalar, Dynamic, 1> storage(stride * n);
+  Map<Mat, 0, OuterStride<>> work(storage.data(), n, n, OuterStride<>(stride));
+  work = a;
+  PartialPivLU<Ref<Mat>> solver(work);
+  const Mat lower = work.template triangularView<UnitLower>();
+  const Mat upper = work.template triangularView<Upper>();
+  const Scalar residual = (solver.permutationP() * a - lower * upper).norm();
+  const Scalar bound = Scalar(64 * n) * NumTraits<Scalar>::epsilon() * a.norm();
+  if (!(residual <= bound)) {
+    state.SkipWithError("LU reconstruction failed");
+    return;
+  }
+  for (auto _ : state) {
+    state.PauseTiming();
+    work = a;
+    state.ResumeTiming();
+    solver.compute(work);
+    benchmark::DoNotOptimize(work.data());
+    benchmark::ClobberMemory();
+  }
+}
+
+BENCHMARK_TEMPLATE(BM_PartialPivLUStride, float, ColMajor)
+    ->ArgsProduct({{32, 128, 500, 512, 768, 1000, 1001, 1024}, {0, 8}});
+BENCHMARK_TEMPLATE(BM_PartialPivLUStride, double, ColMajor)
+    ->ArgsProduct({{32, 128, 500, 512, 768, 1000, 1001, 1024}, {0, 8}});
+BENCHMARK_TEMPLATE(BM_PartialPivLUStride, float, RowMajor)->ArgsProduct({{32, 128, 500, 512, 1000, 1024}, {0, 8}});
+BENCHMARK_TEMPLATE(BM_PartialPivLUStride, double, RowMajor)->ArgsProduct({{32, 128, 500, 512, 1000, 1024}, {0, 8}});
diff --git a/test/householder.cpp b/test/householder.cpp
index e3f9f6f..6c7280b 100644
--- a/test/householder.cpp
+++ b/test/householder.cpp
@@ -783,6 +783,222 @@
   }
 }
 
+template <typename Scalar, int Size, int StorageOrder>
+void householder_short_strided() {
+  using Mat = Matrix<Scalar, Dynamic, Dynamic, StorageOrder>;
+  using RealScalar = typename NumTraits<Scalar>::Real;
+  using Vector = Matrix<Scalar, Size, 1>;
+  using Essential = Matrix<Scalar, Size - 1, 1>;
+  const Vector vector = Vector::Random();
+  Essential essential;
+  Scalar tau;
+  RealScalar beta;
+  vector.makeHouseholder(essential, tau, beta);
+  Vector v;
+  v << Scalar(1), essential;
+  const Matrix<Scalar, Size, Size> h = Matrix<Scalar, Size, Size>::Identity() - tau * v * v.adjoint();
+
+  for (Index cols : {0, 1, 2, 17, 32, 33, 500, 512}) {
+    for (Index innerStride : {1, 2}) {
+      for (Index outerStride : {2048, 2055}) {
+        Matrix<Scalar, Dynamic, 1> storage =
+            Matrix<Scalar, Dynamic, 1>::Random(outerStride * (StorageOrder == ColMajor ? cols + 2 : Size + 2));
+        Map<Mat, 0, Stride<Dynamic, Dynamic>> mapped(storage.data(), Size + 2, cols + 2,
+                                                     Stride<Dynamic, Dynamic>(outerStride, innerStride));
+        const auto originalStorage = storage.eval();
+        const Mat original = mapped;
+        auto block = mapped.block(1, 1, Size, cols);
+        const Mat expected = h * original.block(1, 1, Size, cols);
+        Matrix<Scalar, Dynamic, 1> workspace(cols + 2);
+        workspace.setConstant(Scalar(7));
+        // Fixed essential length selects the fused column-major path even with dynamic block dimensions.
+        block.applyHouseholderOnTheLeft(essential, tau, workspace.data() + 1);
+        const RealScalar bound = RealScalar(16 * Size) * NumTraits<RealScalar>::epsilon() * original.norm();
+        VERIFY((block - expected).norm() <= bound);
+        VERIFY_IS_EQUAL(workspace[0], Scalar(7));
+        VERIFY_IS_EQUAL(workspace[cols + 1], Scalar(7));
+        Mat expectedFull = original;
+        expectedFull.block(1, 1, Size, cols) = block;
+        const auto resultStorage = storage.eval();
+        storage = originalStorage;
+        mapped = expectedFull;
+        VERIFY_IS_EQUAL(storage, resultStorage);
+
+        block.applyHouseholderOnTheLeft(essential, Scalar(0), workspace.data() + 1);
+        VERIFY_IS_EQUAL(storage, resultStorage);
+
+        storage = originalStorage;
+        const Matrix<Scalar, Dynamic, 1> dynamicEssential = essential;
+        block.applyHouseholderOnTheLeft(dynamicEssential, tau, workspace.data() + 1);
+        VERIFY((block - expected).norm() <= bound);
+
+        if (cols > 0) {
+          storage = originalStorage;
+          const Scalar aliasedTau = block.coeff(0, 0);
+          const Matrix<Scalar, Size, Size> aliasedH =
+              Matrix<Scalar, Size, Size>::Identity() - aliasedTau * v * v.adjoint();
+          const Mat aliasedExpected = aliasedH * original.block(1, 1, Size, cols);
+          block.applyHouseholderOnTheLeft(essential, block.coeffRef(0, 0), workspace.data() + 1);
+          VERIFY((block - aliasedExpected).norm() <= bound);
+        }
+      }
+    }
+  }
+}
+
+// Reference for H * a with v = [1; essential] and the operation order of applyHouseholderOnTheLeft, which a
+// noncommutative scalar distinguishes: w_j = tau * (a_0j + sum_i conj(e_i) a_(i+1)j), a_0j -= w_j, a_(i+1)j -= e_i w_j.
+template <typename MatrixType, typename EssentialType>
+MatrixType householder_left_reference(const MatrixType& a, const EssentialType& essential,
+                                      const typename MatrixType::Scalar& tau) {
+  using Scalar = typename MatrixType::Scalar;
+  MatrixType result = a;
+  for (Index j = 0; j < a.cols(); ++j) {
+    Scalar dot = a(0, j);
+    for (Index i = 0; i < essential.size(); ++i) dot += numext::conj(essential(i)) * a(i + 1, j);
+    const Scalar w = tau * dot;
+    result(0, j) -= w;
+    for (Index i = 0; i < essential.size(); ++i) result(i + 1, j) -= essential(i) * w;
+  }
+  return result;
+}
+
+// Whether applyHouseholderOnTheLeft on Derived with this essential type selects the fused column loop.
+template <typename Derived, typename EssentialPart>
+constexpr bool householder_fused() {
+  return std::is_same<internal::householder_apply_left_impl<Derived, EssentialPart>,
+                      internal::householder_apply_left_impl<Derived, EssentialPart, true>>::value;
+}
+
+template <bool ExpectFused, typename MatrixType, typename EssentialType>
+void verify_householder_left(MatrixType a, const EssentialType& essential, const typename MatrixType::Scalar& tau) {
+  using RealScalar = typename NumTraits<typename MatrixType::Scalar>::Real;
+  STATIC_CHECK((householder_fused<MatrixType, EssentialType>() == ExpectFused));
+  const MatrixType expected = householder_left_reference(a, essential.eval(), tau);
+  // Applying a reflector is backward stable; both evaluations carry O(rows * eps * |a|) rounding.
+  const RealScalar bound = RealScalar(16 * a.rows()) * NumTraits<RealScalar>::epsilon() * a.norm();
+  Matrix<typename MatrixType::Scalar, Dynamic, 1> workspace(a.cols());
+  a.applyHouseholderOnTheLeft(essential, tau, workspace.data());
+  VERIFY((a - expected).norm() <= bound);
+}
+
+// The essential part may be any column-vector expression, not only a plain vector: a runtime one-column block
+// (no compile-time vector shape), the adjoint of a row tail as CompleteOrthogonalDecomposition passes, or a
+// fixed-size segment of either. Each must compile in C++14, where no branch of the dispatch is discarded, and
+// only the fixed-size ones on column-major storage select the fused loop.
+template <typename Scalar, int StorageOrder>
+void householder_essential_expressions() {
+  using Mat = Matrix<Scalar, Dynamic, Dynamic, StorageOrder>;
+  using RowMajorMat = Matrix<Scalar, Dynamic, Dynamic, RowMajor>;
+  using VectorType = Matrix<Scalar, Dynamic, 1>;
+  constexpr bool kFusedIfFixed = StorageOrder == ColMajor;
+  const Index cols = internal::random<Index>(1, 20);
+  for (Index rows : {Index(2), Index(3), internal::random<Index>(4, 20)}) {
+    const Mat a = Mat::Random(rows, cols);
+    Scalar tau;
+    typename NumTraits<Scalar>::Real beta;
+    VectorType essential(rows - 1);
+    VectorType::Random(rows).makeHouseholder(essential, tau, beta);
+
+    Mat column(rows + 1, 3);
+    column.col(1).tail(rows - 1) = essential;
+    verify_householder_left<false>(a, column.block(2, 1, rows - 1, 1), tau);
+
+    Mat row(3, rows + 1);
+    row.row(1).tail(rows - 1) = essential.adjoint();
+    verify_householder_left<false>(a, row.row(1).tail(rows - 1).adjoint(), tau);
+
+    RowMajorMat rowMajorColumn(rows + 1, 3);
+    rowMajorColumn.col(1).tail(rows - 1) = essential;
+    verify_householder_left<false>(a, rowMajorColumn.col(1).tail(rows - 1), tau);
+
+    if (rows == 2) {
+      verify_householder_left<kFusedIfFixed>(a, row.row(1).template segment<1>(2).adjoint(), tau);
+      verify_householder_left<kFusedIfFixed>(a, rowMajorColumn.col(1).template tail<1>(), tau);
+    }
+    if (rows == 3) {
+      verify_householder_left<kFusedIfFixed>(a, row.row(1).template segment<2>(2).adjoint(), tau);
+      verify_householder_left<kFusedIfFixed>(a, rowMajorColumn.col(1).template tail<2>(), tau);
+    }
+  }
+}
+
+namespace noncommutative_scalar {
+
+// Quaternions: multiplication is associative and conjugation reverses products, but ab != ba in general.
+struct Quaternion {
+  double r, i, j, k;
+  Quaternion(double real = 0, double x = 0, double y = 0, double z = 0) : r(real), i(x), j(y), k(z) {}
+  Quaternion operator+(const Quaternion& b) const { return Quaternion(r + b.r, i + b.i, j + b.j, k + b.k); }
+  Quaternion operator-(const Quaternion& b) const { return Quaternion(r - b.r, i - b.i, j - b.j, k - b.k); }
+  Quaternion operator-() const { return Quaternion(-r, -i, -j, -k); }
+  Quaternion operator*(const Quaternion& b) const {
+    return Quaternion(r * b.r - i * b.i - j * b.j - k * b.k, r * b.i + i * b.r + j * b.k - k * b.j,
+                      r * b.j - i * b.k + j * b.r + k * b.i, r * b.k + i * b.j - j * b.i + k * b.r);
+  }
+  Quaternion& operator+=(const Quaternion& b) { return *this = *this + b; }
+  Quaternion& operator-=(const Quaternion& b) { return *this = *this - b; }
+  Quaternion& operator*=(const Quaternion& b) { return *this = *this * b; }
+  bool operator==(const Quaternion& b) const { return r == b.r && i == b.i && j == b.j && k == b.k; }
+  bool operator!=(const Quaternion& b) const { return !(*this == b); }
+};
+
+inline Quaternion conj(const Quaternion& a) { return Quaternion(a.r, -a.i, -a.j, -a.k); }
+inline double real(const Quaternion& a) { return a.r; }
+inline double imag(const Quaternion& a) { return a.i; }
+
+}  // namespace noncommutative_scalar
+
+namespace Eigen {
+template <>
+struct NumTraits<noncommutative_scalar::Quaternion> : GenericNumTraits<noncommutative_scalar::Quaternion> {
+  using Real = double;
+  using Literal = double;
+  static constexpr bool IsComplex = true;
+  static constexpr bool RequireInitialization = true;
+};
+}  // namespace Eigen
+
+// The fused loop must multiply in the order the general update's expressions spell out, essential * tmp with
+// tmp = tau * (essential.adjoint() * column + top), which only a noncommutative scalar observes. Small integer
+// components keep every operation exact. The general path is not a reference here: its product kernels reorder
+// operands (the gemv selector transposes essential.adjoint() * bottom, the column-major outer-product selector
+// forms tmp_j * essential).
+template <int Size>
+void householder_noncommutative_scalar() {
+  using Scalar = noncommutative_scalar::Quaternion;
+  using Mat = Matrix<Scalar, Dynamic, Dynamic>;
+  using EssentialType = Matrix<Scalar, Size - 1, 1>;
+  STATIC_CHECK((householder_fused<Mat, EssentialType>()));
+  const Scalar i(0, 1), j(0, 0, 1), k(0, 0, 0, 1);
+  {
+    // With essential = [i, 0...]^T, tau = 1 and input [j, 0...]^T, the second entry becomes -i * j = -k.
+    Mat a = Mat::Constant(Size, 1, Scalar(0));
+    a(0, 0) = j;
+    EssentialType essential = EssentialType::Constant(Scalar(0));
+    essential(0) = i;
+    Scalar workspace[1];
+    a.applyHouseholderOnTheLeft(essential, Scalar(1), workspace);
+    VERIFY(a(0, 0) == Scalar(0));
+    VERIFY(a(1, 0) == -k);
+  }
+  const Index cols = internal::random<Index>(1, 7);
+  Mat a(Size, cols);
+  for (Index c = 0; c < cols; ++c)
+    for (Index r = 0; r < Size; ++r)
+      a(r, c) = Scalar(internal::random<int>(-2, 2), internal::random<int>(-2, 2), internal::random<int>(-2, 2),
+                       internal::random<int>(-2, 2));
+  EssentialType essential;
+  for (Index r = 0; r + 1 < Size; ++r)
+    essential(r) = Scalar(internal::random<int>(-2, 2), internal::random<int>(-2, 2), internal::random<int>(-2, 2),
+                          internal::random<int>(-2, 2));
+  const Scalar tau(internal::random<int>(-2, 2), internal::random<int>(-2, 2), internal::random<int>(-2, 2), 1);
+  const Mat expected = householder_left_reference(a, essential, tau);
+  Matrix<Scalar, Dynamic, 1> workspace(cols);
+  a.applyHouseholderOnTheLeft(essential, tau, workspace.data());
+  VERIFY(a == expected);
+}
+
 EIGEN_DECLARE_TEST(householder) {
   for (int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST_1(householder(Matrix<double, 2, 2>()));
@@ -807,4 +1023,23 @@
   CALL_SUBTEST_11(householder_blocked_right_regression<std::complex<double>>());
   CALL_SUBTEST_12(householder_small_tail());
   CALL_SUBTEST_13(householder_large_components());
+  CALL_SUBTEST_14((householder_short_strided<float, 2, ColMajor>()));
+  CALL_SUBTEST_14((householder_short_strided<float, 3, ColMajor>()));
+  CALL_SUBTEST_15((householder_short_strided<double, 2, ColMajor>()));
+  CALL_SUBTEST_15((householder_short_strided<double, 3, ColMajor>()));
+  CALL_SUBTEST_16((householder_short_strided<std::complex<float>, 2, ColMajor>()));
+  CALL_SUBTEST_16((householder_short_strided<std::complex<float>, 3, ColMajor>()));
+  CALL_SUBTEST_17((householder_short_strided<std::complex<double>, 2, ColMajor>()));
+  CALL_SUBTEST_17((householder_short_strided<std::complex<double>, 3, ColMajor>()));
+  CALL_SUBTEST_18((householder_short_strided<double, 2, RowMajor>()));
+  CALL_SUBTEST_18((householder_short_strided<double, 3, RowMajor>()));
+  CALL_SUBTEST_19((householder_short_strided<std::complex<double>, 3, RowMajor>()));
+  for (int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST_20((householder_essential_expressions<double, ColMajor>()));
+    CALL_SUBTEST_20((householder_essential_expressions<double, RowMajor>()));
+    CALL_SUBTEST_21((householder_essential_expressions<std::complex<double>, ColMajor>()));
+    CALL_SUBTEST_21((householder_essential_expressions<std::complex<double>, RowMajor>()));
+    CALL_SUBTEST_22(householder_noncommutative_scalar<2>());
+    CALL_SUBTEST_22(householder_noncommutative_scalar<3>());
+  }
 }
diff --git a/test/lu.cpp b/test/lu.cpp
index b5e17f0..72131ac 100644
--- a/test/lu.cpp
+++ b/test/lu.cpp
@@ -411,6 +411,36 @@
   }
 }
 
+template <typename Scalar, int StorageOrder>
+void lu_strided_pivots() {
+  using Mat = Matrix<Scalar, Dynamic, Dynamic, StorageOrder>;
+  using RealScalar = typename NumTraits<Scalar>::Real;
+  for (Index n : {17, 32, 33, 127, 128, 129, 256, 512}) {
+    Mat a = Mat::Random(n, n);
+    // A cyclic permutation forces overlapping row swaps across panel boundaries.
+    const Index shift = n / 2 + 1;
+    a.bottomLeftCorner(n - shift, n - shift).diagonal().array() += RealScalar(2 * n);
+    a.topRightCorner(shift, shift).diagonal().array() += RealScalar(2 * n);
+    for (Index padding : {0, 7}) {
+      const Index stride = n + padding;
+      Matrix<Scalar, Dynamic, 1> storage = Matrix<Scalar, Dynamic, 1>::Constant(stride * n + 2, Scalar(7));
+      Map<Mat, 0, OuterStride<>> work(storage.data() + 1, n, n, OuterStride<>(stride));
+      work = a;
+      PartialPivLU<Ref<Mat>> lu(work);
+      VERIFY_IS_EQUAL(lu.matrixLU().data(), work.data());
+      const RealScalar bound = RealScalar(32 * n) * NumTraits<RealScalar>::epsilon() * a.norm();
+      const Mat lower = work.template triangularView<UnitLower>();
+      const Mat upper = work.template triangularView<Upper>();
+      VERIFY((lu.permutationP() * a - lower * upper).norm() <= bound);
+      VERIFY_IS_EQUAL(lu.permutationP().indices()[0], n - shift);
+      VERIFY_IS_EQUAL(storage[0], Scalar(7));
+      VERIFY_IS_EQUAL(storage[storage.size() - 1], Scalar(7));
+      for (Index j = 0; j < n; ++j)
+        for (Index i = n; i < stride; ++i) VERIFY_IS_EQUAL(storage[1 + j * stride + i], Scalar(7));
+    }
+  }
+}
+
 EIGEN_DECLARE_TEST(lu) {
   for (int i = 0; i < g_repeat; i++) {
     CALL_SUBTEST_1(lu_non_invertible<Matrix3f>());
@@ -460,6 +490,13 @@
     CALL_SUBTEST_9(test_2889());
   }
 
+  CALL_SUBTEST_10((lu_strided_pivots<float, ColMajor>()));
+  CALL_SUBTEST_11((lu_strided_pivots<double, ColMajor>()));
+  CALL_SUBTEST_12((lu_strided_pivots<std::complex<float>, ColMajor>()));
+  CALL_SUBTEST_13((lu_strided_pivots<std::complex<double>, ColMajor>()));
+  CALL_SUBTEST_14((lu_strided_pivots<double, RowMajor>()));
+  CALL_SUBTEST_15((lu_strided_pivots<std::complex<double>, RowMajor>()));
+
   // Blocking and vectorization boundary tests (deterministic, outside g_repeat).
   CALL_SUBTEST_3(lu_blocking_boundary<float>());
   CALL_SUBTEST_4(lu_blocking_boundary<double>());