diff --git a/onnxruntime/core/providers/coreml/builders/helper.cc b/onnxruntime/core/providers/coreml/builders/helper.cc
index b8ebbd05a2a20..e1f148fa93e23 100644
--- a/onnxruntime/core/providers/coreml/builders/helper.cc
+++ b/onnxruntime/core/providers/coreml/builders/helper.cc
@@ -50,8 +50,8 @@ bool IsNodeSupported(const Node& node, const OpBuilderInputParams& input_params,
   }
 }
 
-bool IsInputSupported(const Node& node, const NodeArg& input,
-                      const OpBuilderInputParams& input_params, const logging::Logger& logger) {
+bool IsInputSupported(const Node& node, const NodeArg& input, const OpBuilderInputParams& input_params,
+                      const logging::Logger& logger, bool allow_empty_input) {
   if (!input.Exists()) {
     // optional input that is not provided
     return true;
@@ -84,16 +84,10 @@ bool IsInputSupported(const Node& node, const NodeArg& input,
       return false;
     }
 
-    if (dim == 0) {
-      if (node.OpType() == "Resize" && &input == node.InputDefs()[1]) {
-        // one special case. Resize 'roi' input was originally a required input but is rarely used.
-        // ROI is not supported in the CoreML implementation so we will ignore the value, but is often added
-        // (at least in the unit tests) as an initializer with shape {0}.
-      } else {
-        LOGS(logger, WARNING) << "CoreML does not support shapes with dimension values of 0. Input:" << input_name
-                              << ", shape: " << Shape2String(shape);
-        return false;
-      }
+    if (dim == 0 && !allow_empty_input) {
+      LOGS(logger, WARNING) << "CoreML does not support shapes with dimension values of 0. Input:" << input_name
+                            << ", shape: " << Shape2String(shape);
+      return false;
     }
   }
 
diff --git a/onnxruntime/core/providers/coreml/builders/helper.h b/onnxruntime/core/providers/coreml/builders/helper.h
index 300de2dedd122..0acaa0dd8a4a3 100644
--- a/onnxruntime/core/providers/coreml/builders/helper.h
+++ b/onnxruntime/core/providers/coreml/builders/helper.h
@@ -30,7 +30,8 @@ OpBuilderInputParams MakeOpBuilderParams(const GraphViewer& graph_viewer,
 const IOpBuilder* GetOpBuilder(const Node& node);
 
 bool IsInputSupported(const Node& node, const NodeArg& node_arg, const OpBuilderInputParams& input_params,
-                      const logging::Logger& logger);
+                      const logging::Logger& logger,
+                      bool allow_empty_input = false);
 
 bool IsNodeSupported(const Node& node, const OpBuilderInputParams& input_params, const logging::Logger& logger);
 
diff --git a/onnxruntime/core/providers/coreml/builders/impl/base_op_builder.cc b/onnxruntime/core/providers/coreml/builders/impl/base_op_builder.cc
index 83a572f4b60fa..2cae85a0a1c8d 100644
--- a/onnxruntime/core/providers/coreml/builders/impl/base_op_builder.cc
+++ b/onnxruntime/core/providers/coreml/builders/impl/base_op_builder.cc
@@ -74,7 +74,7 @@ bool BaseOpBuilder::IsOpSupported(const Node& node, const OpBuilderInputParams&
 bool BaseOpBuilder::HasSupportedInputs(const Node& node, const OpBuilderInputParams& input_params,
                                        const logging::Logger& logger) const {
   for (const auto* input : node.InputDefs()) {
-    if (!IsInputSupported(node, *input, input_params, logger)) {
+    if (!IsInputSupported(node, *input, input_params, logger, allow_empty_tensor_as_input_)) {
       return false;
     }
   }
diff --git a/onnxruntime/core/providers/coreml/builders/impl/base_op_builder.h b/onnxruntime/core/providers/coreml/builders/impl/base_op_builder.h
index 4a23640d0f34c..071008520fbdc 100644
--- a/onnxruntime/core/providers/coreml/builders/impl/base_op_builder.h
+++ b/onnxruntime/core/providers/coreml/builders/impl/base_op_builder.h
@@ -28,6 +28,10 @@ class BaseOpBuilder : public IOpBuilder {
   void AddInitializersToSkip(ModelBuilder& /*model_builder*/, const Node& /*node*/) const override {}
 
  protected:
+  explicit BaseOpBuilder(bool allow_empty_tensor_as_input = false)
+      : allow_empty_tensor_as_input_(allow_empty_tensor_as_input) {
+  }
+
   // currently we only support float
   static bool IsInputFloat(const Node& node, size_t idx, const OpBuilderInputParams& input_params,
                            const logging::Logger& logger);
@@ -50,6 +54,8 @@ class BaseOpBuilder : public IOpBuilder {
 
   virtual Status AddToModelBuilderImpl(ModelBuilder& model_builder, const Node& node,
                                        const logging::Logger& logger) const = 0;
+
+  const bool allow_empty_tensor_as_input_;  // some operators can handle ignoring an empty tensor as input
 };
 
 }  // namespace coreml
diff --git a/onnxruntime/core/providers/coreml/builders/impl/resize_op_builder.cc b/onnxruntime/core/providers/coreml/builders/impl/resize_op_builder.cc
index 3400f09b4056f..65b5c17f2c6a6 100644
--- a/onnxruntime/core/providers/coreml/builders/impl/resize_op_builder.cc
+++ b/onnxruntime/core/providers/coreml/builders/impl/resize_op_builder.cc
@@ -1,13 +1,15 @@
 // Copyright (c) Microsoft Corporation. All rights reserved.
 // Licensed under the MIT License.
 
-#include <math.h>
+#include <cmath>
 
 #include "core/framework/tensorprotoutils.h"
 #include "core/optimizer/initializer.h"
 #include "core/providers/common.h"
+#include "core/providers/utils.h"
 #include "core/providers/coreml/builders/helper.h"
 #include "core/providers/coreml/builders/impl/base_op_builder.h"
+#include "core/providers/coreml/builders/impl/builder_utils.h"
 #include "core/providers/coreml/builders/model_builder.h"
 #include "core/providers/coreml/builders/op_builder_factory.h"
 #include "core/providers/coreml/shape_utils.h"
@@ -18,6 +20,11 @@ namespace onnxruntime {
 namespace coreml {
 
 class ResizeOpBuilder : public BaseOpBuilder {
+ public:
+  // allow roi and scales potentially being empty inputs that are ignored during processing
+  ResizeOpBuilder() : BaseOpBuilder(/*allow empty inputs*/ true) {}
+
+ private:
   void AddInitializersToSkip(ModelBuilder& model_builder, const Node& node) const override;
 
   Status AddToModelBuilderImpl(ModelBuilder& model_builder, const Node& node,
@@ -29,196 +36,382 @@ class ResizeOpBuilder : public BaseOpBuilder {
   // Resize opset 10- is very different than Resize opset 11+, with many key attributes missing
   // We only support Resize opset 11+ here
   int GetMinSupportedOpSet(const Node& /* node */) const override { return 11; }
+
+  bool SupportsMLProgram() const override { return true; }
 };
 
 namespace {
-bool GetResizeScales(const InitializedTensorSet& initializers,
-                     const Node& node, std::vector<float>& scales,
-                     const logging::Logger&) {
+std::vector<int64_t> GetAxes(const NodeAttrHelper& helper, size_t input_rank) {
+  auto axes = helper.Get("axes", std::vector<int64_t>{});
+  if (axes.empty()) {
+    axes.resize(input_rank);
+    std::iota(axes.begin(), axes.end(), 0);
+  } else {
+    for (auto& value : axes) {
+      if (value < 0) {
+        value = HandleNegativeAxis(value, input_rank);
+      }
+    }
+  }
+
+  return axes;
+}
+
+bool GetValidatedResizeScales(const GraphViewer& graph_viewer,
+                              const Node& node,
+                              const std::vector<int64_t>& input_shape,
+                              const std::vector<int64_t>& axes,
+                              std::vector<float>& scales,
+                              const logging::Logger& logger) {
   const auto& input_defs = node.InputDefs();
-  if (input_defs.size() < 3)
+  int64_t input_rank = input_shape.size();
+
+  if (input_shape[input_rank - 2] == -1 || input_shape[input_rank - 1] == -1) {
+    LOGS(logger, VERBOSE) << "Resize with 'scales' requires the H and W dimensions to have fixed values";
     return false;
+  }
 
-  const auto& scales_tensor = *initializers.at(input_defs[2]->Name());
-  if (scales_tensor.dims_size() != 1 || scales_tensor.dims()[0] != 4)
+  const auto* scales_tensor = graph_viewer.GetConstantInitializer(input_defs[2]->Name());
+  if (!scales_tensor) {
+    LOGS(logger, VERBOSE) << "Resize 'scales' input must be a constant initializer";
     return false;
-  Initializer unpacked_tensor(scales_tensor);
+  }
+
+  Initializer unpacked_tensor(*scales_tensor);
   auto scales_data = unpacked_tensor.DataAsSpan<float>();
-  scales = std::vector<float>{scales_data.begin(), scales_data.end()};
+  scales.assign(scales_data.begin(), scales_data.end());
+
+  for (size_t idx = 0, end = axes.size(); idx < end; ++idx) {
+    auto axis = axes[idx];
+    auto scale = scales[idx];
+    if (axis < (input_rank - 2) && scale != 1.0f) {
+      LOGS(logger, VERBOSE) << "Resize only supports resizing the last two axes. Scale of axis " << axis << " is "
+                            << scale;
+      return false;
+    }
+  }
+
   return true;
 }
 
-bool GetResizeOutputSizes(const InitializedTensorSet& initializers,
-                          const Node& node, std::vector<int64_t>& sizes,
-                          const logging::Logger&) {
+bool GetValidatedResizeSizes(const GraphViewer& graph_viewer,
+                             const Node& node,
+                             const std::vector<int64_t>& input_shape,
+                             const std::vector<int64_t>& axes,
+                             std::vector<int64_t>& sizes, const logging::Logger& logger) {
   const auto& input_defs = node.InputDefs();
-  if (input_defs.size() < 4)
-    return false;
+  int64_t input_rank = input_shape.size();
 
-  const auto& sizes_tensor = *initializers.at(input_defs[3]->Name());
-  if (sizes_tensor.dims_size() != 1 || sizes_tensor.dims()[0] != 4)
+  const auto* sizes_tensor = graph_viewer.GetConstantInitializer(input_defs[3]->Name());
+  if (!sizes_tensor) {
+    LOGS(logger, VERBOSE) << "Resize 'sizes' input must be a constant initializer";
     return false;
-  Initializer unpacked_tensor(sizes_tensor);
+  }
+
+  Initializer unpacked_tensor(*sizes_tensor);
   auto sizes_data = unpacked_tensor.DataAsSpan<int64_t>();
-  sizes = std::vector<int64_t>(sizes_data.begin(), sizes_data.end());
+  sizes.assign(sizes_data.begin(), sizes_data.end());
+
+  for (size_t idx = 0, end = axes.size(); idx < end; ++idx) {
+    auto axis = axes[idx];
+    auto cur_size = input_shape[idx];
+    auto new_size = sizes[idx];
+    if (axis < (input_rank - 2) && cur_size != new_size) {
+      LOGS(logger, VERBOSE) << "Resize only supports resizing the last two axes. Input rank: " << input_rank
+                            << " Change to size of axis " << axis << " from " << cur_size << " to " << new_size;
+      return false;
+    }
+  }
+
   return true;
 }
 }  // namespace
 
 void ResizeOpBuilder::AddInitializersToSkip(ModelBuilder& model_builder, const Node& node) const {
-  // We don't really use ROI here, so add it to skipped list if it's an initializer tensor
-  model_builder.AddInitializerToSkip(node.InputDefs()[1]->Name());  // ROI
-  model_builder.AddInputToSkip(node.InputDefs()[1]->Name());        // ROI
-
-  // We will still add scales to the skipped list even sizes are present
-  // since there is no use of it, we will not process it later
-  model_builder.AddInitializerToSkip(node.InputDefs()[2]->Name());  // scales
-  model_builder.AddInputToSkip(node.InputDefs()[2]->Name());        // scales
-
-  if (node.InputDefs().size() > 3) {
-    model_builder.AddInitializerToSkip(node.InputDefs()[3]->Name());  // sizes
-    model_builder.AddInputToSkip(node.InputDefs()[3]->Name());        // sizes
+  const auto& input_defs = node.InputDefs();
+
+  // In Resize-11 both roi and scales were required even if you were using sizes.
+  // https://github.com/onnx/onnx/blob/main/docs/Changelog.md#Resize-11
+  // From Resize-13 on they're all optional.
+  //
+  // We don't support roi so would never take a node with meaningful roi input. The roi input can however be provided
+  // and is ignored unless coordinate_transformation_mode is set to 'tf_crop_and_resize'.
+  // e.g. our unit tests tend to always provide an empty tensor as roi input instead of as a missing optional input.
+  // Due to this we always call AddInputToSkip on the roi input.
+  //
+  // We require the sizes or scales input to be a constant initializers to take the node (i.e. they won't be an input
+  // to the CoreML model for the partition, so calling AddInputToSkip isn't relevant).
+  // Individual values from scales and sizes are added directly to the layer, so we won't use the initializer.
+  //
+  // That leaves an edge case for Resize-11 where scales could have been provided as an empty input tensor but
+  // we're using a constant initializer for sizes. In this case AddInputToSkip needs to be called for the scales input.
+
+  model_builder.AddInitializerToSkip(input_defs[1]->Name());  // roi
+  model_builder.AddInputToSkip(input_defs[1]->Name());
+
+  if (input_defs[2]->Exists()) {
+    model_builder.AddInitializerToSkip(input_defs[2]->Name());  // scales
+  }
+
+  if (input_defs.size() > 3 && input_defs[3]->Exists()) {
+    model_builder.AddInitializerToSkip(input_defs[3]->Name());  // sizes
+
+    if (node.SinceVersion() < 13) {
+      model_builder.AddInputToSkip(input_defs[2]->Name());  // skip the unused scales input
+    }
   }
 }
 
-Status ResizeOpBuilder::AddToModelBuilderImpl(ModelBuilder& model_builder,
-                                              const Node& node,
+Status ResizeOpBuilder::AddToModelBuilderImpl(ModelBuilder& model_builder, const Node& node,
                                               const logging::Logger& logger) const {
-  std::unique_ptr<COREML_SPEC::NeuralNetworkLayer> layer = model_builder.CreateNNLayer(node);
+  const auto input_defs = node.InputDefs();
+  const auto output_defs = node.OutputDefs();
+  const auto& graph_viewer = model_builder.GetGraphViewer();
+
+  std::vector<int64_t> input_shape;
+  ORT_RETURN_IF_NOT(GetShape(*input_defs[0], input_shape, logger), "Error getting input shape");
+  size_t input_rank = input_shape.size();
+
+  // we know we have either a scales or sizes input so this is safe.
+  // check for sizes first. this handles Resize-11 where scales was a required input but sizes were used if provided.
+  bool using_sizes = input_defs.size() >= 4 && input_defs[3]->Exists();
+  bool using_scales = !using_sizes;
 
-  auto* coreml_upsample = layer->mutable_upsample();
   NodeAttrHelper helper(node);
-  const auto mode = helper.Get("mode", "nearest");
-  if (mode == "linear") {
-    coreml_upsample->set_mode(COREML_SPEC::UpsampleLayerParams_InterpolationMode_BILINEAR);
-  } else {  // we already checked the mode must be NN or Bilinear in IsOpSupportedImpl
-    coreml_upsample->set_mode(COREML_SPEC::UpsampleLayerParams_InterpolationMode_NN);
+  const auto& mode = helper.Get("mode", "nearest");
+  bool is_nearest = mode == "nearest";
+  bool is_linear = !is_nearest;
+
+  auto axes = GetAxes(helper, input_rank);
+  std::vector<float> output_scales;
+  std::vector<int64_t> output_sizes;
+  size_t num_scales = 0;
+  size_t num_sizes = 0;
+
+  if (using_scales) {
+    ORT_RETURN_IF_NOT(GetValidatedResizeScales(graph_viewer, node, input_shape, axes, output_scales, logger),
+                      "Error getting validated scales");
+    num_scales = output_scales.size();
+
+    // special case linear downsample.
+    // the CoreML implementation seems to be flaky and gives different outputs on different OS versions.
+    // use bilinear_resize instead. we check in IsOpSupportedImpl that the downsample input is evenly
+    // divisible by the output size so there's no rounding involved.
+    if (is_linear && (output_scales[num_scales - 1] < 1.f || output_scales[num_scales - 2] < 1.f)) {
+      using_scales = false;
+      using_sizes = true;
+      num_sizes = num_scales;
+      output_sizes = input_shape;
+      // only the last two dims have their size changed
+      output_sizes[input_rank - 2] = static_cast<int64_t>(input_shape[input_rank - 2] * output_scales[num_scales - 2]);
+      output_sizes[input_rank - 1] = static_cast<int64_t>(input_shape[input_rank - 1] * output_scales[num_scales - 1]);
+    }
+  } else {
+    ORT_RETURN_IF_NOT(GetValidatedResizeSizes(graph_viewer, node, input_shape, axes, output_sizes, logger),
+                      "Error getting validated sizes");
+    num_sizes = output_sizes.size();
   }
 
-  const auto& input_defs = node.InputDefs();
-  const auto& initializers(model_builder.GetInitializerTensors());
-
-  if (input_defs.size() >= 3 && input_defs[2]->Exists()) {  // use scales
-    std::vector<float> scales;
-    ORT_RETURN_IF_NOT(GetResizeScales(initializers, node, scales, logger), "Error getting resize scales");
-    coreml_upsample->add_scalingfactor(static_cast<int64_t>(scales[2]));
-    coreml_upsample->add_scalingfactor(static_cast<int64_t>(scales[3]));
-  } else {  // we already checked number of inputs in IsOpSupportedImpl
-    std::vector<int64_t> input_shape;
-    ORT_RETURN_IF_NOT(GetStaticShape(*input_defs[0], input_shape, logger), "Error getting input shape");
-    std::vector<int64_t> output_sizes;
-    ORT_RETURN_IF_NOT(GetResizeOutputSizes(initializers, node, output_sizes, logger),
-                      "Error getting resize output_sizes");
-    coreml_upsample->add_scalingfactor(static_cast<int64_t>(output_sizes[2] / input_shape[2]));
-    coreml_upsample->add_scalingfactor(static_cast<int64_t>(output_sizes[3] / input_shape[3]));
-  }
+#if defined(COREML_ENABLE_MLPROGRAM)
+  if (model_builder.CreateMLProgram()) {
+    using namespace CoreML::Specification::MILSpec;  // NOLINT
+
+    std::string_view coreml_op_type;
+    if (using_scales) {
+      // https://apple.github.io/coremltools/source/coremltools.converters.mil.mil.ops.defs.html#coremltools.converters.mil.mil.ops.defs.iOS15.image_resizing.upsample_bilinear
+      // https://apple.github.io/coremltools/source/coremltools.converters.mil.mil.ops.defs.html#coremltools.converters.mil.mil.ops.defs.iOS15.image_resizing.upsample_nearest_neighbor
+      coreml_op_type = is_linear ? "upsample_bilinear" : "upsample_nearest_neighbor";
+    } else {
+      // https://apple.github.io/coremltools/source/coremltools.converters.mil.mil.ops.defs.html#coremltools.converters.mil.mil.ops.defs.iOS15.image_resizing.resize_bilinear
+      // https://apple.github.io/coremltools/source/coremltools.converters.mil.mil.ops.defs.html#coremltools.converters.mil.mil.ops.defs.iOS15.image_resizing.resize_nearest_neighbor
+      coreml_op_type = is_linear ? "resize_bilinear" : "resize_nearest_neighbor";
+    }
+
+    std::unique_ptr<Operation> op = model_builder.CreateOperation(node, coreml_op_type);
+    AddOperationInput(*op, "x", input_defs[0]->Name());
+
+    std::string coord_trans_mode = helper.Get("coordinate_transformation_mode", "half_pixel");
+
+    if (using_scales) {
+      float scale_height = output_scales[num_scales - 2];
+      float scale_width = output_scales[num_scales - 1];
+      AddOperationInput(*op, "scale_factor_height",
+                        model_builder.AddScalarConstant(coreml_op_type, "scale_factor_height", scale_height));
+      AddOperationInput(*op, "scale_factor_width",
+                        model_builder.AddScalarConstant(coreml_op_type, "scale_factor_width", scale_width));
+
+      if (is_linear) {
+        // we only allow these coord modes in the 'is supported' check,
+        //   - half_pixel or pytorch_half_pixel with output size > 1 -> align_corners = false
+        //   - align_corners -> align_corners = true
+        bool align_corners = coord_trans_mode == "align_corners";
+        AddOperationInput(*op, "align_corners",
+                          model_builder.AddScalarConstant(coreml_op_type, "align_corners", align_corners));
+      }
+    } else {
+      assert(using_sizes);
+      int64_t target_height = output_sizes[num_sizes - 2];
+      int64_t target_width = output_sizes[num_sizes - 1];
+
+      AddOperationInput(*op, "target_size_height",
+                        model_builder.AddScalarConstant(coreml_op_type, "target_size_height", target_height));
+      AddOperationInput(*op, "target_size_width",
+                        model_builder.AddScalarConstant(coreml_op_type, "target_size_width", target_width));
+
+      if (is_linear) {
+        // we only allow these coord modes in the 'is supported' check,
+        //   - half_pixel or pytorch_half_pixel with output size > 1 -> UNALIGN_CORNERS
+        //   - align_corners -> STRICT_ALIGN_CORNERS
+        //   - asymmetric -> DEFAULT
+        std::string sampling_mode_value;
+        if (coord_trans_mode == "asymmetric") {
+          sampling_mode_value = "DEFAULT";
+        } else if (coord_trans_mode == "align_corners") {
+          sampling_mode_value = "STRICT_ALIGN_CORNERS";
+        } else {
+          sampling_mode_value = "UNALIGN_CORNERS";
+        }
+
+        AddOperationInput(*op, "sampling_mode",
+                          model_builder.AddScalarConstant(coreml_op_type, "sampling_mode", sampling_mode_value));
+      }
+    }
 
-  *layer->mutable_input()->Add() = input_defs[0]->Name();
-  *layer->mutable_output()->Add() = node.OutputDefs()[0]->Name();
+    AddOperationOutput(*op, *output_defs[0]);
+    model_builder.AddOperation(std::move(op));
+  } else  // NOLINT
+#endif
+  {
+    std::unique_ptr<COREML_SPEC::NeuralNetworkLayer> layer = model_builder.CreateNNLayer(node);
+
+    auto* coreml_upsample = layer->mutable_upsample();
+
+    // we already checked the mode must be NN or Bilinear in IsOpSupportedImpl
+    if (is_linear) {
+      coreml_upsample->set_mode(COREML_SPEC::UpsampleLayerParams_InterpolationMode_BILINEAR);
+    } else {
+      coreml_upsample->set_mode(COREML_SPEC::UpsampleLayerParams_InterpolationMode_NN);
+    }
+
+    if (using_scales) {
+      coreml_upsample->add_scalingfactor(static_cast<int64_t>(output_scales[num_scales - 2]));
+      coreml_upsample->add_scalingfactor(static_cast<int64_t>(output_scales[num_scales - 1]));
+    } else {
+      auto scale_height = output_sizes[num_sizes - 2] / input_shape[input_rank - 2];
+      auto scale_width = output_sizes[num_sizes - 1] / input_shape[input_rank - 1];
+      coreml_upsample->add_scalingfactor(static_cast<int64_t>(scale_height));
+      coreml_upsample->add_scalingfactor(static_cast<int64_t>(scale_width));
+    }
+
+    *layer->mutable_input()->Add() = input_defs[0]->Name();
+    *layer->mutable_output()->Add() = output_defs[0]->Name();
+
+    model_builder.AddLayer(std::move(layer));
+  }
 
-  model_builder.AddLayer(std::move(layer));
   return Status::OK();
 }
 
 bool ResizeOpBuilder::IsOpSupportedImpl(const Node& node, const OpBuilderInputParams& input_params,
                                         const logging::Logger& logger) const {
   const auto& input_defs = node.InputDefs();
-  const auto& initializers = input_params.graph_viewer.GetAllInitializedTensors();
 
   std::vector<int64_t> input_shape;
-  if (!GetShape(*input_defs[0], input_shape, logger))
+  if (!GetShape(*input_defs[0], input_shape, logger)) {
+    LOGS(logger, VERBOSE) << "Resize: input shape was not known";
     return false;
+  }
 
-  const auto input_size = input_shape.size();
-  if (input_size != 4) {
-    LOGS(logger, VERBOSE) << "Resize only support 4d shape, input is "
-                          << input_size << "d shape";
+  // as we allow empty shapes in the checks done by BaseOpBuilder::HasSupportedInputs we explicitly check for an empty
+  // an empty input here to be consistent.
+  // this should never happen in a real model though as a dim with value 0 (i.e. no input data) would typically be a
+  // dynamic dimension where a previous step had no output (e.g. Loop of zero interations, NonZero with no matches,
+  // NonMaxSupression with no boxes).
+  if (DoesShapeSpecifyZeroElements(input_shape)) {
+    LOGS(logger, VERBOSE) << "Resize input shape has with dimension values of 0 which is not supported.";
     return false;
   }
 
-  {  // check attributes
-    NodeAttrHelper helper(node);
-    const auto mode = helper.Get("mode", "nearest");
-    bool is_linear_resize = mode == "linear";
-    bool is_nearest_resize = mode == "nearest";
-    if (!is_linear_resize && !is_nearest_resize) {
-      LOGS(logger, VERBOSE) << "Resize unsupported input mode, " << mode;
+  const auto input_rank = input_shape.size();
+  if (input_params.create_mlprogram) {
+    if (input_rank < 3 || input_rank > 5) {
+      LOGS(logger, VERBOSE) << "Resize only supports 3D to 5D input. Got: " << input_rank << "D";
       return false;
     }
-
-    const auto exclude_outside = helper.Get("exclude_outside", 0);
-    if (exclude_outside != 0) {
-      LOGS(logger, VERBOSE) << "Resize does not support exclude_outside for now";
+  } else {
+    if (input_rank != 4) {
+      LOGS(logger, VERBOSE) << "Resize only support 4d shape. Got: " << input_rank << "D";
       return false;
     }
+  }
 
-    const auto coord_trans_mode = helper.Get("coordinate_transformation_mode", "half_pixel");
-    bool using_asymmetric = coord_trans_mode == "asymmetric";
-    if (is_linear_resize) {
-      // TODO, add support of align_corners and half_pixel
-      if (!using_asymmetric) {
-        LOGS(logger, VERBOSE) << "Resize bilinear, unsupported coord_trans_mode, " << coord_trans_mode;
-        return false;
-      }
-    } else {
-      // nearest neighbor resizing
-      // For resize using nearest neighbor, we only support coord_trans_mode == "asymmetric" && nearest_mode == "floor"
-      if (!using_asymmetric) {
-        LOGS(logger, VERBOSE) << "Resize nearest neighbor, unsupported coord_trans_mode, " << coord_trans_mode;
-        return false;
-      }
+  // check attributes
+  NodeAttrHelper helper(node);
 
-      const auto nearest_mode = helper.Get("nearest_mode", "round_prefer_floor");
-      if (nearest_mode != "floor") {
-        LOGS(logger, VERBOSE) << "Resize nearest neighbor, unsupported nearest_mode, " << nearest_mode;
-        return false;
-      }
-    }
+  if (helper.Get("antialias", 0) != 0) {
+    LOGS(logger, VERBOSE) << "Resize does not support antialias";
+    return false;
   }
 
-  {  // scales and sizes (if present) must be initializers
-    if (input_defs.size() < 3) {
-      LOGS(logger, VERBOSE) << "Input scales or sizes of Resize must be known";
-      return false;
-    }
+  const auto& mode = helper.Get("mode", "nearest");
+  bool is_linear = mode == "linear";
+  bool is_nearest = mode == "nearest";
+  if (!is_linear && !is_nearest) {
+    LOGS(logger, VERBOSE) << "Resize unsupported input mode: " << mode;
+    return false;
+  }
 
-    bool using_scales = input_defs.size() >= 3 && input_defs[2]->Exists();
-    // scales
-    if (using_scales && !input_params.graph_viewer.GetConstantInitializer(input_defs[2]->Name())) {
-      LOGS(logger, VERBOSE) << "scales input of Resize must be a constant initializer";
+  if (is_nearest) {
+    const auto nearest_mode = helper.Get("nearest_mode", "round_prefer_floor");
+    if (nearest_mode != "floor") {
+      LOGS(logger, VERBOSE) << "Resize only supports 'floor' nearest_mode. Got: " << nearest_mode;
       return false;
     }
+  }
 
-    // sizes
-    if (!using_scales &&
-        (input_defs.size() < 4 ||
-         !input_defs[3]->Exists() ||
-         !input_params.graph_viewer.GetConstantInitializer(input_defs[3]->Name()))) {
-      LOGS(logger, VERBOSE) << "sizes input of Resize must be a constant initializer";
-      return false;
-    }
+  if (helper.Get("exclude_outside", 0) != 0) {
+    LOGS(logger, VERBOSE) << "Resize does not support 'exclude_outside'";
+    return false;
+  }
 
-    // We want to check if the scales or sizes are not trying to resize on N/C channels here
-    if (using_scales) {
-      std::vector<float> scales;
-      if (!GetResizeScales(initializers, node, scales, logger))
-        return false;
+  const auto keep_aspect_ratio_policy = helper.Get("keep_aspect_ratio_policy", "stretch");
+  if (keep_aspect_ratio_policy != "stretch") {
+    LOGS(logger, VERBOSE) << "Resize only supports keep_aspect_ratio_policy of 'stretch'. Got "
+                          << keep_aspect_ratio_policy;
+    return false;
+  }
 
-      float scale_n = scales[0];
-      float scale_c = scales[1];
-      if (scale_n != 1.0f || scale_c != 1.0f) {
-        LOGS(logger, VERBOSE) << "Scales of N/C channel should be 1"
-                              << "Resize of N/C channels are not supported"
-                              << ", scale_n, " << scale_n << ", scale_c, " << scale_c;
-        return false;
-      }
+  // check for sizes first. this handles Resize-11 where scales was a required input but sizes were used if provided.
+  bool using_sizes = input_defs.size() >= 4 && input_defs[3]->Exists();
+  bool using_scales = !using_sizes && input_defs.size() >= 3 && input_defs[2]->Exists();
 
-      // For now we only support upscale, so the scale_h and scale_w should be an integer >= 1
-      // TODO support ResizeBilinear
-      float scale_h = scales[2];
-      float scale_w = scales[3];
+  if (!using_scales && !using_sizes) {
+    LOGS(logger, VERBOSE) << "Resize requires 'scales' or 'sizes' input";
+    return false;
+  }
+
+  // 'axes' is from opset 18 on and allows scales or sizes to have entries for the subset of axes.
+  // we fill with default values if necessary so that the processing is consistent across all supported opsets.
+  auto axes = GetAxes(helper, input_rank);
+  std::vector<float> output_scales;
+  std::vector<int64_t> output_sizes;
+
+  // make sure scales/sizes are constant initializers, and are only modifying the last two dimensions of the input.
+  if (using_scales) {
+    if (!GetValidatedResizeScales(input_params.graph_viewer, node, input_shape, axes, output_scales, logger)) {
+      return false;
+    }
 
-      // Onnx spec requires scale to be a positive float, so we are not checking that here
+    size_t num_scales = output_scales.size();
+    float scale_h = output_scales[num_scales - 2];
+    float scale_w = output_scales[num_scales - 1];
+
+    // NeuralNetwork supports upsample only with round numbers.
+    //
+    // ML Program results seem to match if round numbers are involved. When downsampling the scaling value should be
+    // 1 / <factor of input size>. e.g. if input size is 8, scaling factor could be 1/8, 1/4 or 1/2.
+    if (scale_h >= 1.f && scale_w >= 1.f) {
+      // upsample (or no-op with both == 1.f that we won't bother special-casing)
       if (roundf(scale_h) != scale_h) {
         LOGS(logger, VERBOSE) << "Resize: scale_h: " << scale_h << " is not a whole number";
         return false;
@@ -228,33 +421,57 @@ bool ResizeOpBuilder::IsOpSupportedImpl(const Node& node, const OpBuilderInputPa
         LOGS(logger, VERBOSE) << "Resize: scale_w: " << scale_w << " is not a whole number";
         return false;
       }
-    } else {
-      // we are using sizes
-      std::vector<int64_t> output_sizes;
-      if (!GetResizeOutputSizes(initializers, node, output_sizes, logger))
-        return false;
-
-      if (!IsStaticShape(input_shape)) {
-        LOGS(logger, VERBOSE) << "Input shape with dynamic dimensions is not supported.";
+    } else if (scale_h <= 1.f && scale_w <= 1.f) {
+      // downsample
+      if (input_params.create_mlprogram) {
+        auto h_in = input_shape[input_rank - 2];
+        auto w_in = input_shape[input_rank - 1];
+
+        if (!utils::IsScalingByAFactorOfN(h_in, scale_h)) {
+          LOGS(logger, VERBOSE) << "Resize: downsampling scale " << scale_h
+                                << " is not a factor of input height: " << h_in;
+          return false;
+        }
+
+        if (!utils::IsScalingByAFactorOfN(w_in, scale_w)) {
+          LOGS(logger, VERBOSE) << "Resize: downsampling scale " << scale_w
+                                << " is not a factor of input width: " << w_in;
+          return false;
+        }
+
+      } else {
+        LOGS(logger, VERBOSE) << "Resize: downsampling is not supported.";
         return false;
       }
+    } else {
+      LOGS(logger, VERBOSE) << "Resize: scale_h: " << scale_h << " and scale_w: " << scale_w
+                            << " must both be >= 1 or <= 1";
+      return false;
+    }
+  } else {
+    assert(using_sizes);
+
+    if (!GetValidatedResizeSizes(input_params.graph_viewer, node, input_shape, axes, output_sizes, logger)) {
+      return false;
+    }
 
-      auto output_size_n = output_sizes[0];
-      auto output_size_c = output_sizes[1];
-      if (output_size_n != input_shape[0] || output_size_c != input_shape[1]) {
-        LOGS(logger, VERBOSE) << "Output sizes of N/C channel should match the input sizes, "
-                              << "Resize of N/C channels are not supported"
-                              << ", input_size_n, " << input_shape[0] << ", output_size_n, " << output_size_n
-                              << ". input_size_c, " << input_shape[1] << ", output_size_c, " << output_size_c;
+    if (input_params.create_mlprogram) {
+      // no additional requirements
+    } else {
+      if (!IsStaticShape(input_shape)) {
+        // need to convert from sizes to scales when creating the NN layer, so the input H and W are required
+        LOGS(logger, VERBOSE) << "Resize input shape with dynamic dimensions is not supported.";
         return false;
       }
 
-      // For now we only support upscale, so the output_size_h and output_size_w should be an integer >= 1
+      // For now we only support upsample, so the output_size_h and output_size_w should be an integer >= 1
       // TODO support ResizeBilinear
-      auto output_size_h = output_sizes[2];
-      auto output_size_w = output_sizes[3];
-      auto input_size_h = input_shape[2];
-      auto input_size_w = input_shape[3];
+      auto input_size_h = input_shape[input_rank - 2];
+      auto input_size_w = input_shape[input_rank - 1];
+
+      auto num_sizes = output_sizes.size();  // could be smaller than input_rank if axes was used
+      auto output_size_h = output_sizes[num_sizes - 2];
+      auto output_size_w = output_sizes[num_sizes - 1];
 
       // Onnx spec requires output sizes to be a positive integer, so we are not checking that here
       if (output_size_h % input_size_h != 0) {
@@ -271,6 +488,92 @@ bool ResizeOpBuilder::IsOpSupportedImpl(const Node& node, const OpBuilderInputPa
     }
   }
 
+  std::string coord_trans_mode = helper.Get("coordinate_transformation_mode", "half_pixel");
+  bool using_asymmetric = coord_trans_mode == "asymmetric";
+
+  if (input_params.create_mlprogram) {
+    if (is_nearest) {
+      // Potential CoreML operators we could map to:
+      //
+      // image_resizing.upsample_nearest_neighbor
+      // - mode: nearest
+      // - coordinate_transformation_mode: asymmetric
+      // - 'scales' input
+      //
+      // image_resizing.resize_nearest_neighbor
+      // - mode: nearest
+      // - coordinate_transformation_mode: asymmetric
+      // - 'sizes' input
+      if (!using_asymmetric) {
+        LOGS(logger, VERBOSE) << "Resize with 'mode' of 'nearest' requires 'coordinate_transformation_mode' of "
+                                 "'asymmetric' . Got: "
+                              << coord_trans_mode;
+        return false;
+      }
+    } else {
+      assert(is_linear);
+      // Potential CoreML operators we could map to:
+      //
+      // image_resizing.upsample_bilinear
+      // - mode: linear
+      // - 'scales' input
+      // - coordinate_transformation_mode
+      //   - half_pixel -> align_corners = false
+      //   - align_corners -> align_corners = true
+      //
+      // image_resizing.resize_bilinear
+      // - mode: linear
+      // - 'sizes' input
+      // - coordinate_transformation_mode -> sampling_mode
+      //   - half_pixel -> UNALIGN_CORNERS
+      //   - align_corners -> STRICT_ALIGN_CORNERS
+      //   - asymmetric -> DEFAULT
+      //
+
+      // if output size != 1, coordinate_transformation_mode of pytorch_half_pixel is the same as half_pixel
+      if (coord_trans_mode == "pytorch_half_pixel") {
+        int64_t h_out{0}, w_out{0};
+        if (using_scales) {
+          size_t num_scales = output_scales.size();
+          h_out = std::llround(input_shape[input_rank - 2] * output_scales[num_scales - 2]);
+          w_out = std::llround(input_shape[input_rank - 1] * output_scales[num_scales - 1]);
+        } else {
+          size_t num_sizes = output_sizes.size();
+          h_out = output_sizes[num_sizes - 2];
+          w_out = output_sizes[num_sizes - 1];
+        }
+
+        if (h_out > 1 && w_out > 1) {
+          coord_trans_mode = "half_pixel";
+        }
+      }
+
+      if (coord_trans_mode == "half_pixel" ||
+          coord_trans_mode == "align_corners" ||
+          (using_sizes && coord_trans_mode == "asymmetric")) {
+        // supported
+
+        // FWIW we could calculate (if shape inferencing didn't already) the output sizes and convert a node with
+        // `scales` and co-ord mode of `asymmetric` to having `sizes` input so it's supported.
+      } else {
+        LOGS(logger, VERBOSE) << "Resize with 'mode' of 'linear' requires 'coordinate_transformation_mode' of "
+                                 "'half_pixel', or 'align_corners', or 'asymmetric' with 'sizes' input. Got: "
+                              << coord_trans_mode;
+
+        return false;
+      }
+    }
+  } else {
+    // NeuralNetwork checks
+    if (!using_asymmetric) {
+      // align_corners and half_pixel could be supported in ResizeBilinear but as NeuralNetwork is deprecated
+      // there's no known value to adding that.
+      LOGS(logger, VERBOSE) << "Resize only supports 'asymmetric' coordinate_transformation_mode. Got: "
+                            << coord_trans_mode;
+      return false;
+    }
+  }
+
   return true;
 }
 
diff --git a/onnxruntime/core/providers/coreml/builders/model_builder.h b/onnxruntime/core/providers/coreml/builders/model_builder.h
index 8f85ab2c09e7c..385588dbfdcb8 100644
--- a/onnxruntime/core/providers/coreml/builders/model_builder.h
+++ b/onnxruntime/core/providers/coreml/builders/model_builder.h
@@ -141,8 +141,17 @@ class ModelBuilder {
   // so we don't do a copy of the original initializer into the model.
   void AddInitializerToSkip(const std::string& tensor_name);
 
-  // There are some input which will not be used, add it to a list which will not
-  // be added to CoreML model, since CoreML does not like input unused
+  /// <summary>
+  /// Skip a non-initializer value, that is not used in the CoreML model, but was an input to a supported node.
+  ///
+  /// This is for a rare edge case where a value is an input to a node but is empty/unused, as the
+  /// CoreML model requires all model inputs to be consumed.
+  /// </summary>
+  /// <remarks>
+  /// The only known use case for this currently is Resize, and that is largely due to how the unit tests are
+  /// setup rather than something you'd expect to see in a real model.
+  /// See ResizeOpBuilder::AddInitializersToSkip for more details.
+  /// </remarks>
   void AddInputToSkip(const std::string& input_name);
 
   const std::string& GetUniqueName(const std::string& base_name);
diff --git a/onnxruntime/core/providers/coreml/coreml_execution_provider.cc b/onnxruntime/core/providers/coreml/coreml_execution_provider.cc
index 0ba715cc7c6d9..a92fef81ac395 100644
--- a/onnxruntime/core/providers/coreml/coreml_execution_provider.cc
+++ b/onnxruntime/core/providers/coreml/coreml_execution_provider.cc
@@ -27,6 +27,7 @@ CoreMLExecutionProvider::CoreMLExecutionProvider(uint32_t coreml_flags)
     : IExecutionProvider{onnxruntime::kCoreMLExecutionProvider},
       coreml_flags_(coreml_flags),
       coreml_version_(coreml::util::CoreMLVersion()) {
+  LOGS_DEFAULT(VERBOSE) << "CoreML version: " << coreml_version_;
   if (coreml_version_ < MINIMUM_COREML_VERSION) {
     LOGS_DEFAULT(ERROR) << "CoreML EP is not supported on this platform.";
   }
diff --git a/onnxruntime/core/providers/nnapi/nnapi_builtin/builders/impl/resize_op_builder.cc b/onnxruntime/core/providers/nnapi/nnapi_builtin/builders/impl/resize_op_builder.cc
index d75b9cc72ff4b..ef27f6c942f44 100644
--- a/onnxruntime/core/providers/nnapi/nnapi_builtin/builders/impl/resize_op_builder.cc
+++ b/onnxruntime/core/providers/nnapi/nnapi_builtin/builders/impl/resize_op_builder.cc
@@ -9,6 +9,7 @@
 #include "core/graph/graph_viewer.h"
 #include "core/optimizer/initializer.h"
 #include "core/providers/common.h"
+#include "core/providers/utils.h"
 #include "core/providers/shared/utils/utils.h"
 #include "core/providers/nnapi/nnapi_builtin/builders/helper.h"
 #include "core/providers/nnapi/nnapi_builtin/builders/model_builder.h"
@@ -251,14 +252,34 @@ bool ResizeOpBuilder::IsOpSupportedImpl(const GraphViewer& graph_viewer, const N
       const Initializer unpacked_tensor(*scales);
       auto scales_data = unpacked_tensor.DataAsSpan<float>();
       input_is_nchw = scales_data[1] == 1.0F;
-      float const scale_n = scales_data[0];
-      float const scale_c = input_is_nchw ? scales_data[1] : scales_data[3];
+      const float scale_n = scales_data[0];
+      const float scale_c = input_is_nchw ? scales_data[1] : scales_data[3];
+      const float scale_h = input_is_nchw ? scales_data[2] : scales_data[1];
+      const float scale_w = input_is_nchw ? scales_data[3] : scales_data[2];
+
       if (scale_n != 1.0f || scale_c != 1.0f) {
         LOGS_DEFAULT(VERBOSE) << "Scales of N/C channel should be 1"
                               << "Resize of N/C channels are not supported"
                               << ", scale_n, " << scale_n << ", scale_c, " << scale_c;
         return false;
       }
+
+      // if downsampling the input size must be evenly divisible by the output size to match the onnx output
+      if (scale_h < 1.0f || scale_w < 1.0f) {
+        // we also require input_shape to be known to check
+        auto h_in = input_is_nchw ? input_shape[2] : input_shape[1];
+        auto w_in = input_is_nchw ? input_shape[3] : input_shape[2];
+        if (h_in == 0 || w_in == 0) {
+          LOGS_DEFAULT(VERBOSE) << "Input H and W must be known to downsample with scales";
+          return false;
+        }
+
+        if (!utils::IsScalingByAFactorOfN(h_in, scale_h) ||
+            !utils::IsScalingByAFactorOfN(w_in, scale_w)) {
+          LOGS_DEFAULT(VERBOSE) << "Input size must be evenly divisible by output size when downsampling";
+          return false;
+        }
+      }
     } else {
       const auto* sizes = graph_viewer.GetConstantInitializer(inputs[3].node_arg.Name());
       if (!sizes) {
diff --git a/onnxruntime/core/providers/utils.cc b/onnxruntime/core/providers/utils.cc
index b2f9d265ca053..747b09e42aa21 100644
--- a/onnxruntime/core/providers/utils.cc
+++ b/onnxruntime/core/providers/utils.cc
@@ -23,5 +23,21 @@ common::Status OutputOptionalWithoutDataHelper(const ONNX_NAMESPACE::TypeProto&
   return Status::OK();
 }
 #endif
+
+bool IsScalingByAFactorOfN(int64_t n, float scale) {
+  bool is_factor = false;
+  if (scale > 0.f && scale < 1.f) {
+    const double factor = 1.0 / scale;
+    const double factor_rounded = std::round(factor);
+    constexpr double epsilon = 1.0e-4;  // arbitrarily small enough
+    if (std::abs(factor - factor_rounded) < epsilon) {
+      // result is integer. check if a factor of n
+      const int64_t factor_i = static_cast<int64_t>(factor_rounded);
+      is_factor = n % factor_i == 0;
+    }
+  }
+
+  return is_factor;
+}
 }  // namespace utils
 }  // namespace onnxruntime
diff --git a/onnxruntime/core/providers/utils.h b/onnxruntime/core/providers/utils.h
index 8cafdb8c05cc3..9ea8496a02f85 100644
--- a/onnxruntime/core/providers/utils.h
+++ b/onnxruntime/core/providers/utils.h
@@ -15,5 +15,10 @@ common::Status OutputOptionalWithoutDataHelper(const ONNX_NAMESPACE::TypeProto&
                                                OpKernelContext* context, int output_index);
 #endif
 
+/// <summary>
+/// Check if the reciprocal of 'scale' is a factor of 'n'.
+///   e.g. a scale of 0.5 is 1/2, the reciprocal is 2, and 2 is a factor of any even number.
+/// </summary>
+bool IsScalingByAFactorOfN(int64_t n, float scale);
 }  // namespace utils
 }  // namespace onnxruntime
diff --git a/onnxruntime/core/providers/xnnpack/tensor/resize.cc b/onnxruntime/core/providers/xnnpack/tensor/resize.cc
index 09666c8039402..c752b5f849808 100644
--- a/onnxruntime/core/providers/xnnpack/tensor/resize.cc
+++ b/onnxruntime/core/providers/xnnpack/tensor/resize.cc
@@ -11,6 +11,7 @@
 #include "core/framework/op_kernel.h"
 #include "core/optimizer/initializer.h"
 #include "core/providers/xnnpack/xnnpack_init.h"
+#include "core/providers/utils.h"
 
 namespace onnxruntime {
 namespace xnnpack {
@@ -68,9 +69,27 @@ bool Resize::IsOnnxNodeSupported(const NodeUnit& node_unit,
     InlinedVector<float> scale(4, 1.0F);
     if (scale_tensor) {
       const Initializer scale_val(*scale_tensor, node_unit.ModelPath());
-      if (scale_val.DataAsSpan<float>()[1] != 1.0F) {
+      const auto scales = scale_val.DataAsSpan<float>();
+      if (scales[1] != 1.0F) {
         break;
       }
+
+      // downsampling output seems to require the output size to be a factor of the input to match ONNX
+      if (scales[2] < 1.0f || scales[3] < 1.0f) {
+        // we also require input_shape to be known to check
+        int64_t h_in = x_shape->dim(2).dim_value();
+        int64_t w_in = x_shape->dim(3).dim_value();
+        if (h_in < 0 || w_in < 0) {
+          break;
+        }
+
+        float scale_h = scales[2];
+        float scale_w = scales[3];
+        if (!utils::IsScalingByAFactorOfN(h_in, scale_h) ||
+            !utils::IsScalingByAFactorOfN(w_in, scale_w)) {
+          break;
+        }
+      }
     }
 
     if (size_tensor) {
diff --git a/onnxruntime/test/providers/cpu/tensor/resize_op_test.cc b/onnxruntime/test/providers/cpu/tensor/resize_op_test.cc
index 496f2213e9d32..111520ef03e26 100644
--- a/onnxruntime/test/providers/cpu/tensor/resize_op_test.cc
+++ b/onnxruntime/test/providers/cpu/tensor/resize_op_test.cc
@@ -227,28 +227,33 @@ TEST(ResizeOpTest, NhwcResizeOpLinearDownSampleTest_tf_crop_and_resize_without_e
 }
 
 TEST(ResizeOpTest, ResizeOpLinearDownSampleTest_4DBilinear) {
-  OpTester test("Resize", 13);
-  std::vector<float> roi{};
-  std::vector<float> scales{1.0f, 1.0f, 0.6f, 0.6f};
+  auto run_test = [](bool scales_in_initializer) {
+    OpTester test("Resize", 13);
+    std::vector<float> roi{};
+    std::vector<float> scales{1.0f, 1.0f, 0.6f, 0.6f};
 
-  test.AddAttribute("mode", "linear");
+    test.AddAttribute("mode", "linear");
 
-  constexpr int64_t N = 1, C = 1, H = 2, W = 4;
-  std::vector<float> X = {
-      1.0f, 2.0f, 3.0f, 4.0f,
-      5.0f, 6.0f, 7.0f, 8.0f};
+    constexpr int64_t N = 1, C = 1, H = 2, W = 4;
+    std::vector<float> X = {
+        1.0f, 2.0f, 3.0f, 4.0f,
+        5.0f, 6.0f, 7.0f, 8.0f};
 
-  test.AddInput<float>("X", {N, C, H, W}, X);
-  test.AddInput<float>("roi", {0}, roi);
-  test.AddInput<float>("scales", {4}, scales);
+    test.AddInput<float>("X", {N, C, H, W}, X);
+    test.AddInput<float>("roi", {0}, roi);
+    test.AddInput<float>("scales", {4}, scales, scales_in_initializer);
 
-  std::vector<float> Y = {2.66666651f, 4.3333331f};
+    std::vector<float> Y = {2.66666651f, 4.3333331f};
 
-  test.AddOutput<float>("Y", {N, C, static_cast<int64_t>(H * scales[2]), static_cast<int64_t>(W * scales[3])}, Y);
-  // QNN: result diff
-  // TRT: Segmentation fault in A100
-  std::unordered_set<std::string> excluded_providers({kQnnExecutionProvider});
-  test.Run(OpTester::ExpectResult::kExpectSuccess, "", ExcludeTrtOnA100(excluded_providers));
+    test.AddOutput<float>("Y", {N, C, static_cast<int64_t>(H * scales[2]), static_cast<int64_t>(W * scales[3])}, Y);
+    // QNN: result diff
+    // TRT: Segmentation fault in A100
+    std::unordered_set<std::string> excluded_providers({kQnnExecutionProvider});
+    test.Run(OpTester::ExpectResult::kExpectSuccess, "", ExcludeTrtOnA100(excluded_providers));
+  };
+
+  run_test(false);
+  run_test(true);
 }
 
 TEST(ResizeOpTest, NhwcResizeOpLinearDownSampleTest_4DBilinear) {
@@ -327,13 +332,14 @@ TEST(ResizeOpTest, NhwcResizeOpLinearDownSampleTest_4DBilinear_int8) {
 // Since NNAPI(TFLite) only using the scale calculate using the input/output size
 // For the above test (ResizeOpLinearDownSampleTest_4DBilinear)
 // The output size is [1,1,2,4].*[1,1,0.6,0.6]=[1,1,1,2]
-// NNAPI will recaluclate the scales as the output size divided by input size
+// NNAPI will recalculate the scales as the output size divided by input size
 // scales = [1,1,1,2]./[1,1,2,4] = [1,1,0.5,0.5]
 // See:https://github.com/tensorflow/tensorflow/blob/master/tensorflow/lite/kernels/internal/reference/reference_ops.h
 // So the result of the above example will be different than CPU EP
-// Add the following 2 tests to test with scales valid to NNAPI
+// Add the following 2 tests to test with scales valid to NNAPI.
+// CoreML also doesn't handle a scale that doesn't divide the input size evenly.
 TEST(ResizeOpTest, ResizeOpLinearDownSampleTest_4DBilinear1) {
-  // To test NNAPI EP, we need the sclaes/sizes to be in initializers
+  // To test NNAPI EP, we need the scales/sizes to be in initializers
   auto run_test = [](bool scales_in_initializer) {
     OpTester test("Resize", 13);
     std::vector<float> roi{};
@@ -360,8 +366,38 @@ TEST(ResizeOpTest, ResizeOpLinearDownSampleTest_4DBilinear1) {
   run_test(true);
 }
 
+// Downsize with factor being an odd number (1/3)
+TEST(ResizeOpTest, ResizeOpLinearDownSampleTest_4DBilinear1_OddNumber) {
+  // To test NNAPI EP, we need the scales/sizes to be in initializers
+  auto run_test = [](bool scales_in_initializer) {
+    OpTester test("Resize", 13);
+    std::vector<float> roi{};
+    std::vector<float> scales{1.0f, 1.0f, (1.f / 3), (1.f / 3)};
+
+    test.AddAttribute("mode", "linear");
+
+    constexpr int64_t N = 1, C = 1, H = 3, W = 6;
+    std::vector<float> X = {
+        1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f,
+        7.0f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f,
+        13.0f, 14.0f, 15.0f, 16.0f, 17.0f, 18.0f};
+
+    test.AddInput<float>("X", {N, C, H, W}, X);
+    test.AddInput<float>("roi", {0}, roi);
+    test.AddInput<float>("scales", {4}, scales, scales_in_initializer);
+
+    std::vector<float> Y = {8.f, 11.f};
+
+    test.AddOutput<float>("Y", {N, C, static_cast<int64_t>(H * scales[2]), static_cast<int64_t>(W * scales[3])}, Y);
+    test.Run(OpTester::ExpectResult::kExpectSuccess, "", ExcludeTrtOnA100());
+  };
+
+  run_test(false);
+  run_test(true);
+}
+
 TEST(ResizeOpTest, ResizeOpLinearDownSampleTest_4DBilinear1_WithSizes) {
-  // To test NNAPI EP, we need the sclaes/sizes to be in initializers
+  // To test NNAPI EP, we need the scales/sizes to be in initializers
   auto run_test = [](bool scales_and_sizes_in_initializer) {
     OpTester test("Resize", 13);
     std::vector<float> roi{};
@@ -389,8 +425,32 @@ TEST(ResizeOpTest, ResizeOpLinearDownSampleTest_4DBilinear1_WithSizes) {
   run_test(true);
 }
 
+// test handling for opset 11. scales input is provided but should be ignored in favor of sizes
+TEST(ResizeOpTest, ResizeOpLinearDownSampleTest_4DBilinear1_WithSizesOpset11) {
+  OpTester test("Resize", 11);
+  std::vector<float> roi{};
+  std::vector<float> scales{};
+  constexpr int64_t N = 1, C = 1, H = 2, W = 4;
+  std::vector<int64_t> sizes{N, C, 1, 2};
+  test.AddAttribute("mode", "linear");
+
+  std::vector<float> X = {
+      1.0f, 2.0f, 3.0f, 4.0f,
+      5.0f, 6.0f, 7.0f, 8.0f};
+
+  test.AddInput<float>("X", {N, C, H, W}, X);
+  test.AddInput<float>("roi", {0}, roi);
+  test.AddInput<float>("scales", {0}, scales);
+  test.AddInput<int64_t>("sizes", {4}, sizes, true);  // add as initializer so CoreML EP can take
+
+  std::vector<float> Y = {3.5f, 5.5f};
+
+  test.AddOutput<float>("Y", sizes, Y);
+  test.Run();
+}
+
 TEST(ResizeOpTest, ResizeOpLinearDownSampleTest_4DBilinear_align_corners) {
-  // To test NNAPI EP, we need the sclaes/sizes to be in initializers
+  // To test NNAPI EP, we need the scales/sizes to be in initializers
   auto run_test = [](bool scales_in_initializer) {
     OpTester test("Resize", 13);
     std::vector<float> roi{};
@@ -416,15 +476,51 @@ TEST(ResizeOpTest, ResizeOpLinearDownSampleTest_4DBilinear_align_corners) {
 
   run_test(false);
 
-#ifdef USE_NNAPI
-  // NNAPI will need the scales as an initializer
+#if defined(USE_NNAPI) || defined(USE_COREML)
+  // NNAPI and CoreML need the scales as an initializer
+  // Also tensor RT EP will fail if scales is an initializer but will pass if it is not
+  run_test(true);
+#endif
+}
+
+TEST(ResizeOpTest, ResizeOpLinearDownSampleTest_4DBilinear_align_corners_sizes) {
+  // To test NNAPI EP, we need the scales/sizes to be in initializers
+  auto run_test = [](bool scales_in_initializer) {
+    OpTester test("Resize", 13);
+    std::vector<float> roi{};
+    std::vector<float> scales{};
+    std::vector<int64_t> sizes{1, 1, 1, 2};
+
+    test.AddAttribute("mode", "linear");
+    test.AddAttribute("coordinate_transformation_mode", "align_corners");
+
+    constexpr int64_t N = 1, C = 1, H = 2, W = 4;
+    std::vector<float> X = {
+        1.0f, 2.0f, 3.0f, 4.0f,
+        5.0f, 6.0f, 7.0f, 8.0f};
+
+    test.AddInput<float>("X", {N, C, H, W}, X);
+    test.AddInput<float>("roi", {0}, roi);
+    test.AddInput<float>("", {0}, scales);
+    test.AddInput<int64_t>("sizes", {4}, sizes, scales_in_initializer);
+
+    std::vector<float> Y = {1.0f, 4.0f};
+
+    test.AddOutput<float>("Y", {N, C, 1, 2}, Y);
+    test.Run(OpTester::ExpectResult::kExpectSuccess, "", ExcludeTrtOnA100());
+  };
+
+  run_test(false);
+
+#if defined(USE_NNAPI) || defined(USE_COREML)
+  // NNAPI and CoreML will need the scales as an initializer
   // Also tensor RT EP will fail if scales is an initializer but will pass if it is not
   run_test(true);
 #endif
 }
 
 TEST(ResizeOpTest, NhwcResizeOpLinearDownSampleTest_4DBilinear_align_corners_uint8) {
-  // To test NNAPI EP, we need the sclaes/sizes to be in initializers
+  // To test NNAPI EP, we need the scales/sizes to be in initializers
   auto run_test = [](bool scales_in_initializer) {
     OpTester test("Resize", 13);
     std::vector<float> roi{};
@@ -456,7 +552,7 @@ TEST(ResizeOpTest, NhwcResizeOpLinearDownSampleTest_4DBilinear_align_corners_uin
 }
 
 TEST(ResizeOpTest, NhwcResizeOpLinearDownSampleTest_4DBilinear_align_corners_int8) {
-  // To test NNAPI EP, we need the sclaes/sizes to be in initializers
+  // To test NNAPI EP, we need the scales/sizes to be in initializers
   auto run_test = [](bool scales_in_initializer) {
     OpTester test("Resize", 13);
     std::vector<float> roi{};
@@ -622,7 +718,7 @@ TEST(ResizeOpTest, ResizeOpLinearUpSampleTest_4DBilinear_asymmetric_scales) {
 }
 
 TEST(ResizeOpTest, NhwcResizeOpLinearUpSampleTest_4DBilinear_asymmetric_uint8) {
-  // To test NNAPI EP, we need the sclaes/sizes to be in initializers
+  // To test NNAPI EP, we need the scales/sizes to be in initializers
   auto run_test = [](bool scales_in_initializer) {
     OpTester test("Resize", 13);
     std::vector<float> roi{};
@@ -668,7 +764,7 @@ TEST(ResizeOpTest, NhwcResizeOpLinearUpSampleTest_4DBilinear_asymmetric_uint8) {
 }
 
 TEST(ResizeOpTest, NhwcResizeOpLinearUpSampleTest_4DBilinear_asymmetric_int8) {
-  // To test NNAPI EP, we need the sclaes/sizes to be in initializers
+  // To test NNAPI EP, we need the scales/sizes to be in initializers
   auto run_test = [](bool scales_in_initializer) {
     OpTester test("Resize", 13);
     std::vector<float> roi{};
diff --git a/tools/ci_build/github/apple/coreml_supported_mlprogram_ops.md b/tools/ci_build/github/apple/coreml_supported_mlprogram_ops.md
index 1bbb933f66ba4..3b3790ba06599 100644
--- a/tools/ci_build/github/apple/coreml_supported_mlprogram_ops.md
+++ b/tools/ci_build/github/apple/coreml_supported_mlprogram_ops.md
@@ -17,6 +17,7 @@ Keep in sync with doco generated from /docs/execution-providers/CoreML-Execution
 |ai.onnx:Pow|Only supports cases when both inputs are fp32.|
 |ai.onnx:Relu||
 |ai.onnx:Reshape||
+|ai.onnx:Resize|See [resize_op_builder.cc](https://github.com/microsoft/onnxruntime/blob/main/onnxruntime/core/providers/coreml/builders/impl/resize_op_builder.cc) implementation. There are too many permutations to describe the valid combinations.|
 |ai.onnx:Sub||
 |ai.onnx:Sigmoid||
 |ai:onnx:Tanh||