diff --git a/onnxruntime/core/optimizer/graph_transformer_utils.cc b/onnxruntime/core/optimizer/graph_transformer_utils.cc
index bfbc6fa7aaca0..77315ce4a7041 100644
--- a/onnxruntime/core/optimizer/graph_transformer_utils.cc
+++ b/onnxruntime/core/optimizer/graph_transformer_utils.cc
@@ -57,6 +57,7 @@
 #include "core/optimizer/layer_norm_fusion.h"
 #include "core/optimizer/matmul_activation_fusion.h"
 #include "core/optimizer/matmul_add_fusion.h"
+#include "core/optimizer/multi_head_attention_fusion.h"
 #include "core/optimizer/group_query_attention_pre_norm_fusion.h"
 #include "core/optimizer/matmul_nbits_qkv_fusion.h"
 #include "core/optimizer/matmul_nbits_mlp_fusion.h"
@@ -209,6 +210,26 @@ std::unique_ptr<RuleBasedGraphTransformer> GenerateRuleBasedGraphTransformer(
   return rule_transformer;
 }
 
+#if !defined(DISABLE_CONTRIB_OPS)
+// Whether a registered WebGPU EP has a kernel for the com.microsoft opset-1 operator `op_type`. Fusions that emit
+// WebGPU-only contrib ops check this so they never produce a node the EP cannot run.
+static bool WebGpuEpHasContribKernel(const ExecutionProviders* execution_providers, std::string_view op_type,
+                                     const logging::Logger& logger) {
+  if (execution_providers == nullptr) {
+    return false;
+  }
+  const auto* webgpu_ep = execution_providers->Get(onnxruntime::kWebGpuExecutionProvider);
+  if (webgpu_ep == nullptr) {
+    return false;
+  }
+  auto registry = webgpu_ep->GetKernelRegistry();
+  return registry != nullptr &&
+         registry->TryFindKernel(onnxruntime::kWebGpuExecutionProvider, op_type, kMSDomain, 1,
+                                 KernelRegistry::TypeConstraintMap{}, logger, nullptr)
+             .IsOK();
+}
+#endif  // !defined(DISABLE_CONTRIB_OPS)
+
 InlinedVector<std::unique_ptr<GraphTransformer>> GenerateTransformers(
     TransformerLevel level,
     const SessionOptions& session_options,
@@ -416,6 +437,10 @@ InlinedVector<std::unique_ptr<GraphTransformer>> GenerateTransformers(
       transformers.emplace_back(std::make_unique<LayerNormFusion>(cpu_acl_cuda_dml_webgpu_eps, level));
       transformers.emplace_back(std::make_unique<SimplifiedLayerNormFusion>(cpu_cuda_eps));
       transformers.emplace_back(std::make_unique<AttentionFusion>(cpu_acl_cuda_dml_eps));
+      if (WebGpuEpHasContribKernel(execution_providers, "MultiHeadAttention", logger)) {
+        transformers.emplace_back(std::make_unique<MultiHeadAttentionFusion>(
+            InlinedHashSet<std::string_view>{onnxruntime::kWebGpuExecutionProvider}));
+      }
       transformers.emplace_back(std::make_unique<EmbedLayerNormFusion>(cpu_acl_cuda_dml_eps));
       transformers.emplace_back(std::make_unique<GatherSliceToSplitFusion>(cpu_cuda_eps));
       transformers.emplace_back(std::make_unique<GatherToSliceFusion>(cpu_cuda_eps));
@@ -468,30 +493,12 @@ InlinedVector<std::unique_ptr<GraphTransformer>> GenerateTransformers(
       transformers.emplace_back(std::make_unique<GroupQueryAttentionPreNormFusion>(
           InlinedHashSet<std::string_view>{onnxruntime::kCudaExecutionProvider,
                                            onnxruntime::kWebGpuExecutionProvider}));
-      bool has_matmul_nbits_qlora_kernel = false;
-      bool has_matmul_nbits_mlp_kernel = false;
-      bool has_matmul_nbits_qkv_kernel = false;
-      if (execution_providers != nullptr) {
-        const auto* webgpu_ep = execution_providers->Get(onnxruntime::kWebGpuExecutionProvider);
-        if (webgpu_ep != nullptr) {
-          auto registry = webgpu_ep->GetKernelRegistry();
-          if (registry) {
-            auto has_webgpu_kernel = [&](std::string_view op_type) {
-              return registry->TryFindKernel(onnxruntime::kWebGpuExecutionProvider,
-                                             op_type,
-                                             kMSDomain,
-                                             1,
-                                             KernelRegistry::TypeConstraintMap{},
-                                             logger,
-                                             nullptr)
-                  .IsOK();
-            };
-            has_matmul_nbits_qlora_kernel = has_webgpu_kernel("MatMulNBitsQLoRA");
-            has_matmul_nbits_mlp_kernel = has_webgpu_kernel("MatMulNBitsMlp");
-            has_matmul_nbits_qkv_kernel = has_webgpu_kernel("MatMulNBitsQkv");
-          }
-        }
-      }
+      const bool has_matmul_nbits_qlora_kernel =
+          WebGpuEpHasContribKernel(execution_providers, "MatMulNBitsQLoRA", logger);
+      const bool has_matmul_nbits_mlp_kernel =
+          WebGpuEpHasContribKernel(execution_providers, "MatMulNBitsMlp", logger);
+      const bool has_matmul_nbits_qkv_kernel =
+          WebGpuEpHasContribKernel(execution_providers, "MatMulNBitsQkv", logger);
       if (has_matmul_nbits_qlora_kernel) {
         transformers.emplace_back(std::make_unique<MatMulNBitsQLoRAFusion>());
       }
diff --git a/onnxruntime/core/optimizer/multi_head_attention_fusion.cc b/onnxruntime/core/optimizer/multi_head_attention_fusion.cc
new file mode 100644
index 0000000000000..42b63a6ef0fbe
--- /dev/null
+++ b/onnxruntime/core/optimizer/multi_head_attention_fusion.cc
@@ -0,0 +1,372 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+#include "core/optimizer/multi_head_attention_fusion.h"
+
+#include <array>
+#include <cmath>
+#include <optional>
+
+#include "core/framework/tensorprotoutils.h"
+#include "core/graph/graph_utils.h"
+#include "core/optimizer/initializer.h"
+#include "core/optimizer/utils.h"
+
+namespace onnxruntime {
+
+namespace {
+
+constexpr std::initializer_list<ONNX_NAMESPACE::OperatorSetVersion> kTransposeVersions = {1, 13, 21, 23, 24, 25};
+constexpr std::initializer_list<ONNX_NAMESPACE::OperatorSetVersion> kReshapeVersions = {5, 13, 14, 19, 21, 23, 24, 25};
+constexpr std::initializer_list<ONNX_NAMESPACE::OperatorSetVersion> kMatMulVersions = {1, 9, 13};
+constexpr std::initializer_list<ONNX_NAMESPACE::OperatorSetVersion> kBinaryVersions = {7, 13, 14};
+
+bool IsTransposeWithPerm(const Node& node, const std::vector<int64_t>& perm) {
+  return graph_utils::IsSupportedOptypeVersionAndDomain(node, "Transpose", kTransposeVersions) &&
+         optimizer_utils::IsAttributeWithExpectedValues(node, "perm", perm);
+}
+
+bool IsReshapeWithoutAllowZero(const Node& node) {
+  return graph_utils::IsSupportedOptypeVersionAndDomain(node, "Reshape", kReshapeVersions) &&
+         node.InputDefs().size() == 2 &&
+         !optimizer_utils::IsAttributeWithExpectedValue(node, "allowzero", static_cast<int64_t>(1));
+}
+
+bool GetConstantShape(const Graph& graph, const Node& reshape, InlinedVector<int64_t>& shape) {
+  return optimizer_utils::AppendTensorFromInitializer(graph, *reshape.InputDefs()[1], shape, true);
+}
+
+// Reads a scalar float/float16 constant.
+std::optional<float> GetScalarConstant(const Graph& graph, const NodeArg& arg) {
+  if (!optimizer_utils::IsScalar(arg)) {
+    return std::nullopt;
+  }
+  const auto* tensor_proto = graph_utils::GetConstantInitializer(graph, arg.Name());
+  if (tensor_proto == nullptr) {
+    return std::nullopt;
+  }
+  Initializer init{graph, *tensor_proto, graph.ModelPath()};
+  if (init.size() != 1) {
+    return std::nullopt;
+  }
+  switch (tensor_proto->data_type()) {
+    case ONNX_NAMESPACE::TensorProto_DataType_FLOAT:
+      return *init.data<float>();
+    case ONNX_NAMESPACE::TensorProto_DataType_FLOAT16:
+      return static_cast<float>(*init.data<MLFloat16>());
+    default:
+      return std::nullopt;
+  }
+}
+
+// Matches `Reshape[*,*,N,H] -> Transpose(perm)` feeding `consumer` at `input_index`.
+const Node* MatchHeadSplit(const Graph& graph, const Node& consumer, int input_index, const std::vector<int64_t>& perm,
+                           int64_t& num_heads, int64_t& head_size, const Node*& reshape) {
+  const Node* transpose = graph_utils::GetInputNode(consumer, input_index);
+  if (transpose == nullptr || !IsTransposeWithPerm(*transpose, perm) ||
+      !optimizer_utils::CheckOutputEdges(graph, *transpose, 1)) {
+    return nullptr;
+  }
+
+  reshape = graph_utils::GetInputNode(*transpose, 0);
+  if (reshape == nullptr || !IsReshapeWithoutAllowZero(*reshape)) {
+    return nullptr;
+  }
+
+  InlinedVector<int64_t> shape;
+  if (!GetConstantShape(graph, *reshape, shape) || shape.size() != 4 || shape[2] <= 0 || shape[3] <= 0) {
+    return nullptr;
+  }
+
+  num_heads = shape[2];
+  head_size = shape[3];
+  return transpose;
+}
+
+bool HasRankAtMost(const NodeArg& arg, int rank) {
+  const auto* shape = arg.Shape();
+  return shape != nullptr && shape->dim_size() <= rank;
+}
+
+// Walks up from a head-split Reshape through an optional bias Add and a MatMul to the tensor the projection reads.
+// Only projections that keep the leading dims of that tensor are followed: a 2D constant weight and a bias of rank 1
+// or less, since higher-rank constants could broadcast new leading dims. Gemm projections are not followed: the
+// MatMulAddFusion that produces them requires static shapes, whose batch dims are compared directly.
+const NodeArg* GetProjectionSource(const Graph& graph, const Node& head_split_reshape) {
+  const Node* node = graph_utils::GetInputNode(head_split_reshape, 0);
+  if (node != nullptr && node->OpType() == "Add" && node->InputDefs().size() == 2) {
+    const int data_index = graph_utils::IsConstantInitializer(graph, node->InputDefs()[1]->Name(), true) ? 0 : 1;
+    const NodeArg& bias = *node->InputDefs()[1 - data_index];
+    if (!graph_utils::IsConstantInitializer(graph, bias.Name(), true) || !HasRankAtMost(bias, 1)) {
+      return nullptr;
+    }
+    node = graph_utils::GetInputNode(*node, data_index);
+  }
+  if (node == nullptr || !graph_utils::IsSupportedOptypeVersionAndDomain(*node, "MatMul", kMatMulVersions)) {
+    return nullptr;
+  }
+  const NodeArg& weight = *node->InputDefs()[1];
+  if (!graph_utils::IsConstantInitializer(graph, weight.Name(), true) || !weight.Shape() ||
+      weight.Shape()->dim_size() != 2) {
+    return nullptr;
+  }
+  return node->InputDefs()[0];
+}
+
+// Whether two head-split Reshapes produce the same batch dim. Statically comparable inferred dims decide it. When
+// they are unknown (dynamic axes), it still holds if both project the same tensor through leading-dim-preserving
+// projections and keep the same leading reshape dims, as in self-attention.
+bool HaveSameBatchDim(const Graph& graph, const Node& a, const Node& b) {
+  const auto* a_shape = a.OutputDefs()[0]->Shape();
+  const auto* b_shape = b.OutputDefs()[0]->Shape();
+  if (a_shape != nullptr && b_shape != nullptr && a_shape->dim_size() == 4 && b_shape->dim_size() == 4) {
+    const auto& a_dim = a_shape->dim(0);
+    const auto& b_dim = b_shape->dim(0);
+    if (utils::HasDimValue(a_dim) && utils::HasDimValue(b_dim)) {
+      return a_dim.dim_value() == b_dim.dim_value();
+    }
+    if (utils::HasDimParam(a_dim) && utils::HasDimParam(b_dim) && a_dim.dim_param() == b_dim.dim_param()) {
+      return true;
+    }
+  }
+
+  InlinedVector<int64_t> a_dims, b_dims;
+  const NodeArg* a_source = GetProjectionSource(graph, a);
+  return a_source != nullptr && a_source == GetProjectionSource(graph, b) && GetConstantShape(graph, a, a_dims) &&
+         GetConstantShape(graph, b, b_dims) && a_dims[0] == b_dims[0] && a_dims[1] == b_dims[1];
+}
+
+bool TryFuse(Graph& graph, Node& softmax, const InlinedHashSet<std::string_view>& compatible_eps,
+             NodeArg*& bsd_shape_arg, const logging::Logger& logger) {
+  if (!graph_utils::IsSupportedOptypeVersionAndDomain(softmax, "Softmax", {1, 11, 13}) ||
+      !graph_utils::IsSupportedProvider(softmax, compatible_eps) ||
+      !optimizer_utils::CheckOutputEdges(graph, softmax, 1)) {
+    return false;
+  }
+
+  // Before opset 13 Softmax coerces the input to 2D at `axis` (default 1), so only the last axis of the 4D score
+  // tensor matches MultiHeadAttention.
+  const auto* axis_attr = graph_utils::GetNodeAttribute(softmax, "axis");
+  const int64_t axis = axis_attr != nullptr ? axis_attr->i() : (softmax.SinceVersion() >= 13 ? -1 : 1);
+  if (axis != -1 && axis != 3) {
+    return false;
+  }
+
+  float scale = 1.0f;
+  const Node* scale_node = nullptr;
+  const Node* qk_matmul = graph_utils::GetInputNode(softmax, 0);
+  if (qk_matmul != nullptr &&
+      (graph_utils::IsSupportedOptypeVersionAndDomain(*qk_matmul, "Div", kBinaryVersions) ||
+       graph_utils::IsSupportedOptypeVersionAndDomain(*qk_matmul, "Mul", kBinaryVersions))) {
+    scale_node = qk_matmul;
+    const bool is_div = scale_node->OpType() == "Div";
+    if (!optimizer_utils::CheckOutputEdges(graph, *scale_node, 1)) {
+      return false;
+    }
+    // For Mul the scalar may be either operand; for Div it must be the divisor.
+    int scalar_index = 1;
+    std::optional<float> value = GetScalarConstant(graph, *scale_node->InputDefs()[1]);
+    if (!value.has_value() && !is_div) {
+      scalar_index = 0;
+      value = GetScalarConstant(graph, *scale_node->InputDefs()[0]);
+    }
+    if (!value.has_value() || *value == 0.0f) {
+      return false;
+    }
+    scale = is_div ? 1.0f / *value : *value;
+    // Also catches a tiny divisor whose reciprocal overflows.
+    if (!std::isfinite(scale)) {
+      return false;
+    }
+    qk_matmul = graph_utils::GetInputNode(*scale_node, 1 - scalar_index);
+  }
+
+  if (qk_matmul == nullptr || !graph_utils::IsSupportedOptypeVersionAndDomain(*qk_matmul, "MatMul", kMatMulVersions) ||
+      !optimizer_utils::CheckOutputEdges(graph, *qk_matmul, 1)) {
+    return false;
+  }
+
+  Node* pv_matmul = graph.GetNode(softmax.OutputNodesBegin()->Index());
+  if (pv_matmul == nullptr || !graph_utils::IsSupportedOptypeVersionAndDomain(*pv_matmul, "MatMul", kMatMulVersions) ||
+      pv_matmul->InputDefs()[0] != softmax.OutputDefs()[0] ||
+      !optimizer_utils::CheckOutputEdges(graph, *pv_matmul, 1)) {
+    return false;
+  }
+
+  Node* out_transpose = graph.GetNode(pv_matmul->OutputNodesBegin()->Index());
+  if (out_transpose == nullptr || !IsTransposeWithPerm(*out_transpose, {0, 2, 1, 3}) ||
+      !optimizer_utils::CheckOutputEdges(graph, *out_transpose, 1)) {
+    return false;
+  }
+
+  // The output Reshape is kept and now reads the [B, S, N*H] MultiHeadAttention output instead of the
+  // [B, S, N, H] Transpose output. Both tensors hold the same elements in the same order, so explicit dims and -1
+  // keep their meaning, but a 0 copies the input dim at that index, which only matches for the two leading dims.
+  // The shape therefore has to be a constant without zeros past index 1.
+  const Node* out_reshape = graph.GetNode(out_transpose->OutputNodesBegin()->Index());
+  InlinedVector<int64_t> out_shape;
+  if (out_reshape == nullptr || !IsReshapeWithoutAllowZero(*out_reshape) ||
+      out_reshape->InputDefs()[0] != out_transpose->OutputDefs()[0] ||
+      !GetConstantShape(graph, *out_reshape, out_shape)) {
+    return false;
+  }
+  for (size_t i = 2; i < out_shape.size(); ++i) {
+    if (out_shape[i] == 0) {
+      return false;
+    }
+  }
+
+  int64_t q_heads = 0, q_head_size = 0, k_heads = 0, k_head_size = 0, v_heads = 0, v_head_size = 0;
+  const Node *q_reshape = nullptr, *k_reshape = nullptr, *v_reshape = nullptr;
+  const Node* q_transpose = MatchHeadSplit(graph, *qk_matmul, 0, {0, 2, 1, 3}, q_heads, q_head_size, q_reshape);
+  const Node* k_transpose = MatchHeadSplit(graph, *qk_matmul, 1, {0, 2, 3, 1}, k_heads, k_head_size, k_reshape);
+  const Node* v_transpose = MatchHeadSplit(graph, *pv_matmul, 1, {0, 2, 1, 3}, v_heads, v_head_size, v_reshape);
+  if (q_transpose == nullptr || k_transpose == nullptr || v_transpose == nullptr) {
+    return false;
+  }
+  // Unequal head counts also cover multi-query attention, where K/V have one head that MatMul broadcasts.
+  if (q_heads != k_heads || q_heads != v_heads || q_head_size != k_head_size || q_head_size != v_head_size) {
+    LOGS(logger, VERBOSE) << "MultiHeadAttentionFusion[" << softmax.Name() << "]: mismatched Q/K/V head shapes";
+    return false;
+  }
+  // MatMul broadcasts the batch dim, but MultiHeadAttention requires Q, K and V to share it.
+  if (!HaveSameBatchDim(graph, *q_reshape, *k_reshape) || !HaveSameBatchDim(graph, *q_reshape, *v_reshape)) {
+    return false;
+  }
+
+  const auto* q_type = q_reshape->OutputDefs()[0]->TypeAsProto();
+  if (q_type == nullptr || !q_type->has_tensor_type() ||
+      (q_type->tensor_type().elem_type() != ONNX_NAMESPACE::TensorProto_DataType_FLOAT &&
+       q_type->tensor_type().elem_type() != ONNX_NAMESPACE::TensorProto_DataType_FLOAT16)) {
+    return false;
+  }
+
+  // The remaining checks keep the fused node on the WebGPU kernel's dense FlashAttention prefill path. Its fallbacks
+  // are slower than the unfused graph, and the non-flash one does not support graph capture. They mirror kernel logic
+  // in contrib_ops/webgpu/bert/flash_attention.cc and contrib_ops/cpu/bert/multihead_attention_helper.h, so keep them
+  // in sync with the functions named below.
+  // multihead_attention_helper::CheckInputs rejects more heads than maxComputeInvocationsPerWorkgroup, whose
+  // guaranteed minimum is 256.
+  if (q_heads > 256) {
+    return false;
+  }
+  // CanApplyFlashAttention requires a head size that is a multiple of 8 on every adapter.
+  if (q_head_size % 8 != 0) {
+    return false;
+  }
+  // CanApplyDensePrefill needs K/V tiles of 16 keys to fit in the 16 KB of workgroup storage every adapter
+  // guarantees (DensePrefillFitsWorkgroupStorage, SelectDensePrefillMaxKStep). That matches the kernel except on
+  // Qualcomm, whose extra output tile can exceed 16 KB on 16 KB adapters; the adapter is
+  // unknown here, and missing dense prefill only costs speed, so the bound is not tightened for that case.
+  const int64_t element_size =
+      q_type->tensor_type().elem_type() == ONNX_NAMESPACE::TensorProto_DataType_FLOAT16 ? 2 : 4;
+  if (2 * 16 * q_head_size * element_size > 16384) {
+    return false;
+  }
+  // CanApplyDensePrefill sends queries shorter than 32 to the split-reduce decode kernel, and they leave only a small
+  // score tensor to avoid.
+  // A symbolic query length cannot be checked here, so it is fused.
+  const auto* q_bnsh_shape = q_transpose->OutputDefs()[0]->Shape();
+  if (q_bnsh_shape != nullptr && q_bnsh_shape->dim_size() == 4 && utils::HasDimValue(q_bnsh_shape->dim(2)) &&
+      q_bnsh_shape->dim(2).dim_value() < 32) {
+    return false;
+  }
+
+  const auto& provider = softmax.GetExecutionProviderType();
+  for (const Node* node : {qk_matmul, scale_node, static_cast<const Node*>(pv_matmul),
+                           static_cast<const Node*>(out_transpose), q_transpose, k_transpose, v_transpose,
+                           q_reshape, k_reshape, v_reshape, out_reshape}) {
+    if (node != nullptr && node->GetExecutionProviderType() != provider) {
+      return false;
+    }
+  }
+
+  if (bsd_shape_arg == nullptr) {
+    ONNX_NAMESPACE::TensorProto shape_proto;
+    shape_proto.set_name(graph.GenerateNodeArgName("mha_fusion_bsd_shape"));
+    shape_proto.set_data_type(ONNX_NAMESPACE::TensorProto_DataType_INT64);
+    shape_proto.add_dims(3);
+    const std::array<int64_t, 3> bsd_shape{0, 0, -1};
+    utils::SetRawDataInTensorProto(shape_proto, bsd_shape.data(), bsd_shape.size() * sizeof(int64_t));
+    bsd_shape_arg = &graph_utils::AddInitializerWithOrtValue(graph, shape_proto);
+  }
+
+  // Only the element type is set; shape inference fills in the shapes when the graph is resolved.
+  ONNX_NAMESPACE::TypeProto tensor_type;
+  tensor_type.mutable_tensor_type()->set_elem_type(q_type->tensor_type().elem_type());
+
+  auto add_to_bsd = [&](const Node& reshape) -> Node& {
+    NodeArg* bsnh = graph.GetNode(reshape.Index())->MutableOutputDefs()[0];
+    NodeArg& bsd = graph.GetOrCreateNodeArg(graph.GenerateNodeArgName(bsnh->Name() + "_bsd"), &tensor_type);
+    Node& node = graph.AddNode(graph.GenerateNodeName(reshape.Name() + "/MhaBsd"), "Reshape",
+                               "Fold heads into the hidden dim for MultiHeadAttention", {bsnh, bsd_shape_arg}, {&bsd});
+    node.SetExecutionProviderType(provider);
+    graph.AddEdge(reshape.Index(), node.Index(), 0, 0);
+    return node;
+  };
+  Node& q_bsd = add_to_bsd(*q_reshape);
+  Node& k_bsd = add_to_bsd(*k_reshape);
+  Node& v_bsd = add_to_bsd(*v_reshape);
+
+  NodeArg& mha_output = graph.GetOrCreateNodeArg(graph.GenerateNodeArgName("mha_fusion_output"), &tensor_type);
+  Node& mha = graph.AddNode(graph.GenerateNodeName(softmax.Name() + "/MultiHeadAttention"), "MultiHeadAttention",
+                            "Fused scaled dot-product attention",
+                            {q_bsd.MutableOutputDefs()[0], k_bsd.MutableOutputDefs()[0], v_bsd.MutableOutputDefs()[0]},
+                            {&mha_output}, nullptr, kMSDomain);
+  mha.AddAttribute("num_heads", q_heads);
+  mha.AddAttribute("scale", scale);
+  mha.SetExecutionProviderType(provider);
+  graph.AddEdge(q_bsd.Index(), mha.Index(), 0, 0);
+  graph.AddEdge(k_bsd.Index(), mha.Index(), 0, 1);
+  graph.AddEdge(v_bsd.Index(), mha.Index(), 0, 2);
+
+  graph_utils::ReplaceDownstreamNodeInput(graph, *out_transpose, 0, mha, 0);
+
+  InlinedVector<NodeIndex> nodes_to_remove{q_transpose->Index(), k_transpose->Index(), v_transpose->Index(),
+                                           qk_matmul->Index(), softmax.Index(), pv_matmul->Index(),
+                                           out_transpose->Index()};
+  if (scale_node != nullptr) {
+    nodes_to_remove.push_back(scale_node->Index());
+  }
+  graph.NotifyNodeReplacement(nodes_to_remove, mha.Index());
+  for (NodeIndex index : nodes_to_remove) {
+    Node* node = graph.GetNode(index);
+    graph_utils::RemoveNodeOutputEdges(graph, *node);
+    graph.RemoveNode(index);
+  }
+
+  LOGS(logger, VERBOSE) << "MultiHeadAttentionFusion: fused attention at " << mha.Name();
+  return true;
+}
+
+}  // namespace
+
+Status MultiHeadAttentionFusion::ApplyImpl(Graph& graph, bool& modified, int graph_level,
+                                           const logging::Logger& logger) const {
+  GraphViewer graph_viewer(graph);
+  const auto& node_topology_list = graph_viewer.GetNodesInTopologicalOrder();
+
+  NodeArg* bsd_shape_arg = nullptr;
+  int fused_count = 0;
+  for (auto node_index : node_topology_list) {
+    Node* node = graph.GetNode(node_index);
+    if (node == nullptr) {
+      continue;
+    }
+
+    ORT_RETURN_IF_ERROR(Recurse(*node, modified, graph_level, logger));
+
+    if (TryFuse(graph, *node, GetCompatibleExecutionProviders(), bsd_shape_arg, logger)) {
+      ++fused_count;
+      modified = true;
+    }
+  }
+
+  if (fused_count > 0) {
+    LOGS(logger, INFO) << "Total fused MultiHeadAttention node count: " << fused_count;
+  }
+
+  return Status::OK();
+}
+
+}  // namespace onnxruntime
diff --git a/onnxruntime/core/optimizer/multi_head_attention_fusion.h b/onnxruntime/core/optimizer/multi_head_attention_fusion.h
new file mode 100644
index 0000000000000..eb4f6466f5665
--- /dev/null
+++ b/onnxruntime/core/optimizer/multi_head_attention_fusion.h
@@ -0,0 +1,47 @@
+// Copyright (c) Microsoft Corporation. All rights reserved.
+// Licensed under the MIT License.
+
+#pragma once
+
+#include "core/optimizer/graph_transformer.h"
+
+namespace onnxruntime {
+
+/**
+@Class MultiHeadAttentionFusion
+
+Fuses the unmasked scaled-dot-product attention emitted by ViT-style exporters (e.g. HuggingFace
+ViT / DINOv2) with separate Q, K and V projections into a com.microsoft MultiHeadAttention node:
+
+    Q -> Reshape[*,*,N,H] -> Transpose(0,2,1,3) --+
+                                                  MatMul -> Div|Mul(scalar)? -> Softmax(-1)
+    K -> Reshape[*,*,N,H] -> Transpose(0,2,3,1) --+                                 |
+                                                                                    MatMul -> Transpose(0,2,1,3)
+    V -> Reshape[*,*,N,H] -> Transpose(0,2,1,3) ------------------------------------+         -> Reshape(out)
+
+is rewritten to
+
+    Q -> Reshape[*,*,N,H] -> Reshape[0,0,-1] --+
+    K -> Reshape[*,*,N,H] -> Reshape[0,0,-1] --+-> MultiHeadAttention(num_heads=N, scale) -> Reshape(out)
+    V -> Reshape[*,*,N,H] -> Reshape[0,0,-1] --+
+
+The added Reshapes only fold the trailing (N, H) dims back to the BSD layout MultiHeadAttention expects, so they are
+free for EPs whose Reshape aliases its input. This keeps the [B, N, S, S] score tensor out of memory, which dominates
+the attention cost for long vision sequences.
+
+Besides the structural match, the fusion requires Q, K and V to share their batch dim, head count and head size, and
+it only fires when the WebGPU MultiHeadAttention kernel would take its dense FlashAttention prefill path: a head size
+that is a multiple of 8 and fits the guaranteed workgroup storage, and a query length that is not statically known to
+be below 32. These cutoffs are specific to the WebGPU kernel, so the transformer is only meant to be registered for
+the WebGPU EP.
+*/
+class MultiHeadAttentionFusion : public GraphTransformer {
+ public:
+  explicit MultiHeadAttentionFusion(const InlinedHashSet<std::string_view>& compatible_execution_providers) noexcept
+      : GraphTransformer("MultiHeadAttentionFusion", compatible_execution_providers) {
+  }
+
+  Status ApplyImpl(Graph& graph, bool& modified, int graph_level, const logging::Logger& logger) const override;
+};
+
+}  // namespace onnxruntime
diff --git a/onnxruntime/test/optimizer/graph_transform_test.cc b/onnxruntime/test/optimizer/graph_transform_test.cc
index ec69411e11b6a..813555fcb4f8f 100644
--- a/onnxruntime/test/optimizer/graph_transform_test.cc
+++ b/onnxruntime/test/optimizer/graph_transform_test.cc
@@ -8,6 +8,7 @@
 #include <cmath>
 #include <functional>
 #include <fstream>
+#include <optional>
 #include <random>
 #include <sstream>
 
@@ -30,6 +31,7 @@
 #include "core/graph/model.h"
 #include "core/mlas/inc/mlas_q4.h"
 #include "core/optimizer/attention_fusion.h"
+#include "core/optimizer/multi_head_attention_fusion.h"
 #include "core/optimizer/bias_dropout_fusion.h"
 #include "core/optimizer/bias_gelu_fusion.h"
 #include "core/optimizer/bias_softmax_fusion.h"
@@ -9716,6 +9718,442 @@ TEST_F(GraphTransformationTests, AttentionFusionMobileClipMhaProjectionRewriteFa
                                         CheckMobileClipAttentionUnfusedMatMulGraph));
 }
 
+#if !defined(DISABLE_CONTRIB_OPS)
+enum class VitAttentionScaleOp {
+  Div,
+  Mul,
+  // MatMul feeds Softmax directly.
+  None,
+};
+
+struct VitAttentionTestConfig {
+  VitAttentionScaleOp scale_op = VitAttentionScaleOp::Div;
+  // Mul only: put the scalar in input 0.
+  bool scale_is_first_input = false;
+  // Overrides the Div/Mul scalar, which defaults to the usual sqrt(head_size) divisor or its reciprocal.
+  std::optional<float> scale_value;
+  // Unset means no Softmax axis attribute.
+  std::optional<int64_t> softmax_axis = 3;
+  std::vector<int64_t> k_perm = {0, 2, 3, 1};
+  int64_t seq_len = 64;
+  int64_t num_heads = 4;
+  int64_t head_size = 16;
+  // Fewer K/V heads than the query heads makes MatMul broadcast them, as in multi-query attention.
+  // 0 means the same as num_heads.
+  int64_t kv_num_heads = 0;
+  // Projects K/V from a second input with batch 1; a larger query batch then makes MatMul broadcast them.
+  bool separate_kv_input = false;
+  int64_t q_batch = 1;
+  // Leaves the batch and sequence dims of the inputs unknown, as in models exported with dynamic axes.
+  bool dynamic_input_dims = false;
+  // A non-zero value gives the Q weight this leading batch dim, which MatMul broadcasts into the Q batch.
+  int64_t q_weight_batch = 0;
+  // Adds bias Adds after the projections: Q uses this shape, K/V use their 1D hidden size.
+  std::vector<int64_t> q_bias_shape;
+  // Shape of the Reshape that merges the heads back; empty means {0, 0, hidden_size}.
+  std::vector<int64_t> merge_heads_shape;
+  bool fp16 = false;
+};
+
+// Unmasked ViT/DINOv2 self-attention with separate Q, K and V projections:
+//   x -> MatMul(Wq/Wk/Wv) -> Reshape[0,0,N,H] -> Transpose -> MatMul(Q, K^T) -> Div|Mul -> Softmax
+//     -> MatMul(V) -> Transpose(0,2,1,3) -> Reshape[0,0,N*H]
+// The fp16 variant uses constant data and is only meant for graph checks.
+static void BuildVitAttentionTestCase(ModelTestBuilder& builder, const VitAttentionTestConfig& config = {}) {
+  const int64_t num_heads = config.num_heads;
+  const int64_t kv_num_heads = config.kv_num_heads == 0 ? num_heads : config.kv_num_heads;
+  const int64_t head_size = config.head_size;
+  const int64_t hidden_size = num_heads * head_size;
+
+  auto make_input = [&](int64_t batch) -> NodeArg* {
+    if (config.dynamic_input_dims) {
+      return builder.MakeInput<float>(std::optional<std::vector<int64_t>>{{-1, -1, hidden_size}});
+    }
+    const std::vector<int64_t> shape{batch, config.seq_len, hidden_size};
+    if (config.fp16) {
+      return builder.MakeInput<MLFloat16>(std::optional<std::vector<int64_t>>{shape});
+    }
+    return builder.MakeInput<float>(shape, -1.0f, 1.0f);
+  };
+  auto make_weight = [&](int64_t output_size, int64_t batch = 0) -> NodeArg* {
+    std::vector<int64_t> shape{hidden_size, output_size};
+    if (batch > 0) {
+      shape.insert(shape.begin(), batch);
+    }
+    if (config.fp16) {
+      return builder.MakeInitializer<MLFloat16>(
+          shape, std::vector<MLFloat16>(static_cast<size_t>(TensorShape(shape).Size()), MLFloat16(0.01f)));
+    }
+    return builder.MakeInitializer<float>(shape, -0.2f, 0.2f);
+  };
+
+  NodeArg* q_input = make_input(config.q_batch);
+  NodeArg* kv_input = config.separate_kv_input ? make_input(1) : q_input;
+  const float sqrt_head_size = std::sqrt(static_cast<float>(head_size));
+  const float scale_value = config.scale_value.value_or(
+      config.scale_op == VitAttentionScaleOp::Div ? sqrt_head_size : 1.0f / sqrt_head_size);
+  NodeArg* scale = config.fp16 ? builder.MakeScalarInitializer<MLFloat16>(MLFloat16(scale_value))
+                               : builder.MakeScalarInitializer<float>(scale_value);
+  NodeArg* merge_heads_shape = builder.Make1DInitializer<int64_t>(
+      config.merge_heads_shape.empty() ? std::vector<int64_t>{0, 0, hidden_size} : config.merge_heads_shape);
+
+  auto add_head_split = [&](NodeArg* x, int64_t heads, const std::vector<int64_t>& perm, int64_t weight_batch,
+                            const std::vector<int64_t>& bias_shape) {
+    auto* weight = make_weight(heads * head_size, weight_batch);
+    auto* split_heads_shape = builder.Make1DInitializer<int64_t>({0, 0, heads, head_size});
+    auto* projected = builder.MakeIntermediate();
+    auto* bsnh = builder.MakeIntermediate();
+    auto* transposed = builder.MakeIntermediate();
+    builder.AddNode("MatMul", {x, weight}, {projected});
+    if (!config.q_bias_shape.empty()) {
+      auto* bias = builder.MakeInitializer<float>(bias_shape, -0.1f, 0.1f);
+      auto* biased = builder.MakeIntermediate();
+      builder.AddNode("Add", {bias, projected}, {biased});
+      projected = biased;
+    }
+    builder.AddNode("Reshape", {projected, split_heads_shape}, {bsnh});
+    builder.AddNode("Transpose", {bsnh}, {transposed}).AddAttribute("perm", perm);
+    return transposed;
+  };
+  const std::vector<int64_t> kv_bias_shape{kv_num_heads * head_size};
+  auto* q = add_head_split(q_input, num_heads, {0, 2, 1, 3}, config.q_weight_batch, config.q_bias_shape);
+  auto* k = add_head_split(kv_input, kv_num_heads, config.k_perm, 0, kv_bias_shape);
+  auto* v = add_head_split(kv_input, kv_num_heads, {0, 2, 1, 3}, 0, kv_bias_shape);
+
+  auto* qk = builder.MakeIntermediate();
+  auto* scaled = builder.MakeIntermediate();
+  auto* probs = builder.MakeIntermediate();
+  auto* context = builder.MakeIntermediate();
+  auto* context_bsnh = builder.MakeIntermediate();
+  auto* output = builder.MakeOutput();
+  builder.AddNode("MatMul", {q, k}, {qk});
+  if (config.scale_op == VitAttentionScaleOp::Div) {
+    builder.AddNode("Div", {qk, scale}, {scaled});
+  } else if (config.scale_op == VitAttentionScaleOp::None) {
+    scaled = qk;
+  } else if (config.scale_is_first_input) {
+    builder.AddNode("Mul", {scale, qk}, {scaled});
+  } else {
+    builder.AddNode("Mul", {qk, scale}, {scaled});
+  }
+  auto& softmax = builder.AddNode("Softmax", {scaled}, {probs});
+  if (config.softmax_axis.has_value()) {
+    softmax.AddAttribute("axis", *config.softmax_axis);
+  }
+  builder.AddNode("MatMul", {probs, v}, {context});
+  builder.AddNode("Transpose", {context}, {context_bsnh}).AddAttribute("perm", std::vector<int64_t>{0, 2, 1, 3});
+  builder.AddNode("Reshape", {context_bsnh, merge_heads_shape}, {output});
+}
+
+static Status AssignVitAttentionEp(Graph& graph, const char* ep) {
+  for (auto& node : graph.Nodes()) {
+    node.SetExecutionProviderType(ep);
+  }
+  return Status::OK();
+}
+
+static Status CheckVitAttentionFusedWithScale(Graph& graph, float expected_scale) {
+  auto op_to_count = CountOpsInGraph(graph);
+  ORT_RETURN_IF_NOT(op_to_count["com.microsoft.MultiHeadAttention"] == 1, "expected one MultiHeadAttention");
+  ORT_RETURN_IF_NOT(op_to_count["Softmax"] == 0, "Softmax was not fused");
+  ORT_RETURN_IF_NOT(op_to_count["Transpose"] == 0, "head Transposes were not removed");
+  ORT_RETURN_IF_NOT(op_to_count["MatMul"] == 3, "only the Q/K/V projections should remain");
+  ORT_RETURN_IF_NOT(op_to_count["Div"] == 0 && op_to_count["Mul"] == 0, "the scale was not folded");
+  // Three head-split Reshapes, the three added BSD Reshapes and the output Reshape.
+  ORT_RETURN_IF_NOT(op_to_count["Reshape"] == 7, "unexpected Reshape count");
+
+  for (const auto& node : graph.Nodes()) {
+    if (node.OpType() == "MultiHeadAttention") {
+      const auto& attrs = node.GetAttributes();
+      ORT_RETURN_IF_NOT(attrs.at("num_heads").i() == 4, "unexpected num_heads");
+      ORT_RETURN_IF_NOT(std::abs(attrs.at("scale").f() - expected_scale) < 1e-6f, "unexpected scale");
+      ORT_RETURN_IF_NOT(node.GetExecutionProviderType() == kWebGpuExecutionProvider, "unexpected EP");
+      // Q, K and V must be the [0, 0, -1] Reshapes of the head-split [B, S, N, H] tensors, not those tensors.
+      for (int i = 0; i < 3; ++i) {
+        const Node* bsd = graph_utils::GetInputNode(node, i);
+        ORT_RETURN_IF_NOT(bsd != nullptr && bsd->OpType() == "Reshape", "MultiHeadAttention input ", i,
+                          " is not a BSD Reshape");
+        InlinedVector<int64_t> bsd_shape;
+        const InlinedVector<int64_t> expected_bsd_shape{0, 0, -1};
+        const bool is_bsd_shape =
+            optimizer_utils::AppendTensorFromInitializer(graph, *bsd->InputDefs()[1], bsd_shape, true) &&
+            bsd_shape == expected_bsd_shape;
+        ORT_RETURN_IF_NOT(is_bsd_shape, "MultiHeadAttention input ", i, " is not reshaped to [0, 0, -1]");
+        const Node* head_split = graph_utils::GetInputNode(*bsd, 0);
+        ORT_RETURN_IF_NOT(head_split != nullptr && head_split->OpType() == "Reshape", "MultiHeadAttention input ", i,
+                          " does not come from a head-split Reshape");
+      }
+    }
+  }
+  return Status::OK();
+}
+
+static Status CheckVitAttentionFused(Graph& graph) {
+  return CheckVitAttentionFusedWithScale(graph, 0.25f);
+}
+
+static Status CheckVitAttentionFusedUnscaled(Graph& graph) {
+  return CheckVitAttentionFusedWithScale(graph, 1.0f);
+}
+
+static Status CheckVitAttentionUnfused(Graph& graph) {
+  auto op_to_count = CountOpsInGraph(graph);
+  ORT_RETURN_IF_NOT(op_to_count["com.microsoft.MultiHeadAttention"] == 0, "unexpected MultiHeadAttention");
+  ORT_RETURN_IF_NOT(op_to_count["Softmax"] == 1, "Softmax should remain");
+  return Status::OK();
+}
+
+static Status AssignVitAttentionWebGpuEp(Graph& graph) {
+  return AssignVitAttentionEp(graph, kWebGpuExecutionProvider);
+}
+
+static Status TestVitAttentionFusion(const VitAttentionTestConfig& config, int opset,
+                                     const std::function<Status(Graph&)>& check, const logging::Logger& logger,
+                                     const std::function<Status(Graph&)>& assign_ep = AssignVitAttentionWebGpuEp) {
+  return TestGraphTransformer(
+      [&config](ModelTestBuilder& builder) { BuildVitAttentionTestCase(builder, config); }, opset, logger,
+      std::make_unique<MultiHeadAttentionFusion>(InlinedHashSet<std::string_view>{kWebGpuExecutionProvider}),
+      TransformerLevel::Level2, 1, assign_ep, check);
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitDivScale) {
+  // Opset 11 matches exported ViT/DINOv2 models, where Softmax needs the explicit last-axis attribute.
+  ASSERT_STATUS_OK(TestVitAttentionFusion({}, 11, CheckVitAttentionFused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitMulScaleDefaultAxis) {
+  // From opset 13 Softmax defaults to the last axis.
+  VitAttentionTestConfig config;
+  config.scale_op = VitAttentionScaleOp::Mul;
+  config.softmax_axis = std::nullopt;
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 13, CheckVitAttentionFused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitMulScaleFirstInput) {
+  VitAttentionTestConfig config;
+  config.scale_op = VitAttentionScaleOp::Mul;
+  config.scale_is_first_input = true;
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 11, CheckVitAttentionFused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitNoScale) {
+  VitAttentionTestConfig config;
+  config.scale_op = VitAttentionScaleOp::None;
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 11, CheckVitAttentionFusedUnscaled, *logger_));
+}
+
+// Fails when ONNX bumps the since_version of an op in the pattern beyond the version lists in
+// onnxruntime/core/optimizer/multi_head_attention_fusion.cc.
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionCurrentOpsetTest) {
+  const int current_opset = GetCurrentOnnxOpset();
+  auto check = [current_opset](Graph& graph) -> Status {
+    Status status = CheckVitAttentionFused(graph);
+    if (!status.IsOK()) {
+      return ORT_MAKE_STATUS(ONNXRUNTIME, FAIL, "MultiHeadAttentionFusion failed at opset ", current_opset, ": ",
+                             status.ErrorMessage(), ". Update the version lists in ",
+                             "onnxruntime/core/optimizer/multi_head_attention_fusion.cc.");
+    }
+    return Status::OK();
+  };
+  // As in GeluFusionCurrentOpsetTest, allow the current opset even while it is still under development in ONNX.
+  ASSERT_STATUS_OK(TestGraphTransformer(
+      [](ModelTestBuilder& builder) { BuildVitAttentionTestCase(builder); }, current_opset, *logger_,
+      std::make_unique<MultiHeadAttentionFusion>(InlinedHashSet<std::string_view>{kWebGpuExecutionProvider}),
+      TransformerLevel::Level2, 1, AssignVitAttentionWebGpuEp, check,
+      ModelOptions{kAllowReleasedOpsetsOnly, /*strict_shape_type_inference*/ false}));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitOverflowingScaleNotFused) {
+  // The reciprocal of a denormal divisor overflows to inf, so there is no finite scale attribute to emit.
+  VitAttentionTestConfig config;
+  config.scale_value = 1e-40f;
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 11, CheckVitAttentionUnfused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitFloat16) {
+  VitAttentionTestConfig config;
+  config.fp16 = true;
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 11, CheckVitAttentionFused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitOpset11DefaultAxisNotFused) {
+  // Before opset 13 the default Softmax axis is 1, which normalizes over heads and keys together.
+  VitAttentionTestConfig config;
+  config.softmax_axis = std::nullopt;
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 11, CheckVitAttentionUnfused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitUntransposedKeyNotFused) {
+  // seq_len == head_size keeps Q x K a valid MatMul although it no longer computes Q x K^T.
+  VitAttentionTestConfig config;
+  config.k_perm = {0, 2, 1, 3};
+  config.seq_len = 32;
+  config.head_size = 32;
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 11, CheckVitAttentionUnfused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitMultiQueryNotFused) {
+  VitAttentionTestConfig config;
+  config.kv_num_heads = 1;
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 11, CheckVitAttentionUnfused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitBroadcastBatchNotFused) {
+  VitAttentionTestConfig config;
+  config.separate_kv_input = true;
+  config.q_batch = 2;
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 11, CheckVitAttentionUnfused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitDynamicDimsSelfAttention) {
+  // The inferred batch dims are unknown, but Q, K and V project the same tensor, so they share them.
+  VitAttentionTestConfig config;
+  config.dynamic_input_dims = true;
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 11, CheckVitAttentionFused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitDynamicDimsCrossAttentionNotFused) {
+  // With unknown dims and a separate K/V input, nothing shows that the batch dims match.
+  VitAttentionTestConfig config;
+  config.dynamic_input_dims = true;
+  config.separate_kv_input = true;
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 11, CheckVitAttentionUnfused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitDynamicDimsWithBiasFused) {
+  // 1D biases keep the leading dims, so the shared projection source still proves matching batches.
+  VitAttentionTestConfig config;
+  config.dynamic_input_dims = true;
+  config.q_bias_shape = {64};
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 11, CheckVitAttentionFused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitSharedSourceBatchedWeightNotFused) {
+  // Q, K and V project the same [1, S, D] input, but the [2, D, D] Q weight broadcasts Q to batch 2. The static
+  // inferred batch dims (2 vs. 1) reject it before the projection-source check is consulted.
+  VitAttentionTestConfig config;
+  config.q_weight_batch = 2;
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 11, CheckVitAttentionUnfused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitSharedSourceBatchedWeightDynamicDimsNotFused) {
+  // As above with unknown inferred batch dims, so only the projection-source check can reject it.
+  VitAttentionTestConfig config;
+  config.dynamic_input_dims = true;
+  config.q_weight_batch = 2;
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 11, CheckVitAttentionUnfused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitSharedSourceBatchedBiasDynamicDimsNotFused) {
+  // A [2, 1, D] Q bias broadcasts Q to batch 2 even though all three project the same input.
+  VitAttentionTestConfig config;
+  config.dynamic_input_dims = true;
+  config.q_bias_shape = {2, 1, 64};
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 11, CheckVitAttentionUnfused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitZeroInMergedHeadDimNotFused) {
+  // The output Reshape would read the MultiHeadAttention output, whose dim 2 is N*H rather than N.
+  VitAttentionTestConfig config;
+  config.merge_heads_shape = {0, 0, 0, -1};
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 11, CheckVitAttentionUnfused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitTooManyHeadsNotFused) {
+  // The WebGPU kernel rejects more heads than maxComputeInvocationsPerWorkgroup, which is only guaranteed to be 256.
+  // The other guards pass: the head size is a multiple of 8 and the query is not short.
+  VitAttentionTestConfig config;
+  config.num_heads = 257;
+  config.head_size = 8;
+  config.seq_len = 32;
+  config.fp16 = true;
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 11, CheckVitAttentionUnfused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitShortQueryNotFused) {
+  // Queries shorter than 32 would run the split-reduce decode kernel.
+  VitAttentionTestConfig config;
+  config.seq_len = 31;
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 11, CheckVitAttentionUnfused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitHeadSizeNotMultipleOf8NotFused) {
+  // Such head sizes can miss the FlashAttention path, whose fallback does not support graph capture.
+  VitAttentionTestConfig config;
+  config.head_size = 12;
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 11, CheckVitAttentionUnfused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitLargeFloatHeadNotFused) {
+  // fp32 K/V tiles for a head size of 256 exceed the guaranteed 16 KB of workgroup storage.
+  VitAttentionTestConfig config;
+  config.head_size = 256;
+  ASSERT_STATUS_OK(TestVitAttentionFusion(config, 11, CheckVitAttentionUnfused, *logger_));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitOtherEpNotFused) {
+  ASSERT_STATUS_OK(TestVitAttentionFusion({}, 11, CheckVitAttentionUnfused, *logger_, [](Graph& graph) {
+    return AssignVitAttentionEp(graph, kCpuExecutionProvider);
+  }));
+}
+
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitMixedEpNotFused) {
+  // One head-split Transpose left on the CPU EP.
+  ASSERT_STATUS_OK(TestVitAttentionFusion({}, 11, CheckVitAttentionUnfused, *logger_, [](Graph& graph) {
+    ORT_RETURN_IF_ERROR(AssignVitAttentionEp(graph, kWebGpuExecutionProvider));
+    for (auto& node : graph.Nodes()) {
+      if (node.OpType() == "Transpose") {
+        node.SetExecutionProviderType(kCpuExecutionProvider);
+        break;
+      }
+    }
+    return Status::OK();
+  }));
+}
+
+// Runs the graph on WebGPU with and without the fusion through the default session optimizers and compares the
+// outputs. A sequence length of at least 32 takes MultiHeadAttention down the dense FlashAttention prefill path; 97
+// also leaves a partial K/V tile. Each session gets its own EP instance.
+TEST_F(GraphTransformationTests, MultiHeadAttentionFusionVitWebGpuNumerics) {
+  if (!DefaultWebGpuExecutionProvider()) {
+    GTEST_SKIP() << "WebGPU execution provider is not available";
+  }
+
+  const std::unordered_map<std::string, int> domain_to_version{{kOnnxDomain, 11}, {kMSDomain, 1}};
+  Model model("VitAttention", false, ModelMetaData(), PathString(), IOnnxRuntimeOpSchemaRegistryList(),
+              domain_to_version, {}, *logger_);
+  ModelTestBuilder builder(model.MainGraph());
+  VitAttentionTestConfig config;
+  config.seq_len = 97;
+  BuildVitAttentionTestCase(builder, config);
+  builder.SetGraphOutputs();
+  ASSERT_STATUS_OK(model.MainGraph().Resolve());
+  std::string model_data;
+  model.ToProto().SerializeToString(&model_data);
+
+  auto run = [&](bool fuse, std::vector<OrtValue>& fetches) {
+    InferenceSessionWrapper session{SessionOptions{}, GetEnvironment()};
+    ASSERT_STATUS_OK(session.RegisterExecutionProvider(DefaultWebGpuExecutionProvider()));
+    std::istringstream model_istream(model_data);
+    ASSERT_STATUS_OK(session.Load(model_istream));
+    if (!fuse) {
+      ASSERT_STATUS_OK(session.FilterEnabledOptimizers({"MultiHeadAttentionFusion"}));
+    }
+    ASSERT_STATUS_OK(session.Initialize());
+    auto op_to_count = CountOpsInGraph(session.GetGraph());
+    EXPECT_EQ(op_to_count["com.microsoft.MultiHeadAttention"], fuse ? 1 : 0);
+    ASSERT_STATUS_OK(session.Run(RunOptions{}, builder.feeds_, builder.output_names_, &fetches));
+  };
+
+  std::vector<OrtValue> unfused, fused;
+  ASSERT_NO_FATAL_FAILURE(run(false, unfused));
+  ASSERT_NO_FATAL_FAILURE(run(true, fused));
+  ASSERT_EQ(fused.size(), 1u);
+  auto [result, message] = CompareOrtValue(fused[0], unfused[0], 1e-4, 1e-4, false);
+  EXPECT_EQ(result, COMPARE_RESULT::SUCCESS) << message;
+}
+#endif  // !defined(DISABLE_CONTRIB_OPS)
+
 // Current-opset regression tests for fusion optimizers.
 // These construct minimal graphs at the current ONNX opset and verify the optimizer fires.
 // If ONNX bumps an op's since_version, the optimizer's version list will miss it and the test fails.