diff --git a/include/xtensor/views/xview.hpp b/include/xtensor/views/xview.hpp
index c647572de..6bc9a44d0 100644
--- a/include/xtensor/views/xview.hpp
+++ b/include/xtensor/views/xview.hpp
@@ -817,25 +817,14 @@ namespace xt
     private:
 
         bool is_newaxis_slice(size_type index) const noexcept;
-        void to_end_impl(layout_type l);
-
-        template <class F>
-        void common_step_forward(size_type dim, F f);
-        template <class F>
-        void common_step_backward(size_type dim, F f);
-
-        template <class F>
-        void common_step_forward(size_type dim, size_type n, F f);
-        template <class F>
-        void common_step_backward(size_type dim, size_type n, F f);
-
-        template <class F>
-        void common_reset(size_type dim, F f, bool backwards);
+        void to_slice_begin();
+        void common_reset(size_type dim, bool backwards);
 
         view_type* p_view;
         substepper_type m_it;
         size_type m_offset;
-        std::array<std::size_t, sizeof...(S)> m_index_keeper;
+        // Position of the stepper in each slice, indexed like the slices.
+        std::array<std::size_t, sizeof...(S)> m_index_keeper{};
     };
 
     // meta-function returning the shape type for an xview
@@ -2010,26 +1999,13 @@ namespace xt
         , m_it(it)
         , m_offset(offset)
     {
-        if (!end)
+        if (end)
         {
-            std::fill(m_index_keeper.begin(), m_index_keeper.end(), 0);
-            auto func = [](const auto& s) noexcept
-            {
-                return xt::value(s, 0);
-            };
-            for (size_type i = 0; i < sizeof...(S); ++i)
-            {
-                if (!is_newaxis_slice(i))
-                {
-                    size_type s = apply<size_type>(i, func, p_view->slices());
-                    size_type index = i - newaxis_count_before<S...>(i);
-                    m_it.step(index, s);
-                }
-            }
+            to_end(l);
         }
         else
         {
-            to_end_impl(l);
+            to_slice_begin();
         }
     }
 
@@ -2042,75 +2018,112 @@ namespace xt
     template <bool is_const, class CT, class... S>
     inline void xview_stepper<is_const, CT, S...>::step(size_type dim)
     {
-        auto func = [this](size_type index, size_type offset)
-        {
-            m_it.step(index, offset);
-        };
-        common_step_forward(dim, func);
+        step(dim, size_type(1));
     }
 
     template <bool is_const, class CT, class... S>
     inline void xview_stepper<is_const, CT, S...>::step_back(size_type dim)
     {
-        auto func = [this](size_type index, size_type offset)
-        {
-            m_it.step_back(index, offset);
-        };
-        common_step_backward(dim, func);
+        step_back(dim, size_type(1));
     }
 
+    // dim is a dimension of the (possibly broadcast) shape being iterated; the
+    // m_offset leading dimensions of that shape are not dimensions of the view.
     template <bool is_const, class CT, class... S>
     inline void xview_stepper<is_const, CT, S...>::step(size_type dim, size_type n)
     {
-        auto func = [this](size_type index, size_type offset)
+        if (dim < m_offset)
         {
-            m_it.step(index, offset);
-        };
-        common_step_forward(dim, n, func);
+            return;
+        }
+        const size_type index = integral_skip<S...>(dim - m_offset);
+        if (index >= sizeof...(S))
+        {
+            m_it.step(index - newaxis_count<S...>(), n);
+        }
+        else if (!is_newaxis_slice(index))
+        {
+            auto func = [this, index, n](const auto& s) noexcept
+            {
+                const size_type st = step_size(s, m_index_keeper[index], n);
+                m_index_keeper[index] += n;
+                return st;
+            };
+            m_it.step(index - newaxis_count_before<S...>(index), apply<size_type>(index, func, p_view->slices()));
+        }
     }
 
     template <bool is_const, class CT, class... S>
     inline void xview_stepper<is_const, CT, S...>::step_back(size_type dim, size_type n)
     {
-        auto func = [this](size_type index, size_type offset)
+        if (dim < m_offset)
         {
-            m_it.step_back(index, offset);
-        };
-        common_step_backward(dim, n, func);
+            return;
+        }
+        const size_type index = integral_skip<S...>(dim - m_offset);
+        if (index >= sizeof...(S))
+        {
+            m_it.step_back(index - newaxis_count<S...>(), n);
+        }
+        else if (!is_newaxis_slice(index))
+        {
+            auto func = [this, index, n](const auto& s) noexcept
+            {
+                m_index_keeper[index] -= n;
+                return step_size(s, m_index_keeper[index], n);
+            };
+            m_it.step_back(
+                index - newaxis_count_before<S...>(index),
+                apply<size_type>(index, func, p_view->slices())
+            );
+        }
     }
 
     template <bool is_const, class CT, class... S>
     inline void xview_stepper<is_const, CT, S...>::reset(size_type dim)
     {
-        auto func = [this](size_type index, size_type offset)
-        {
-            m_it.step_back(index, offset);
-        };
-        common_reset(dim, func, false);
+        common_reset(dim, false);
     }
 
     template <bool is_const, class CT, class... S>
     inline void xview_stepper<is_const, CT, S...>::reset_back(size_type dim)
     {
-        auto func = [this](size_type index, size_type offset)
-        {
-            m_it.step(index, offset);
-        };
-        common_reset(dim, func, true);
+        common_reset(dim, true);
     }
 
     template <bool is_const, class CT, class... S>
     inline void xview_stepper<is_const, CT, S...>::to_begin()
     {
-        std::fill(m_index_keeper.begin(), m_index_keeper.end(), 0);
         m_it.to_begin();
+        to_slice_begin();
     }
 
+    // Leaves the stepper where increment_stepper would leave it after the last
+    // element: one step past it along the innermost dimension of the traversal.
     template <bool is_const, class CT, class... S>
     inline void xview_stepper<is_const, CT, S...>::to_end(layout_type l)
     {
-        m_it.to_end(l);
-        to_end_impl(l);
+        if (l != layout_type::row_major && l != layout_type::column_major)
+        {
+            XTENSOR_THROW(std::runtime_error, "Iteration only allowed in row or column major.");
+        }
+        to_begin();
+        if (p_view->size() == 0)
+        {
+            return;
+        }
+        const auto& shape = p_view->shape();
+        const size_type dimension = shape.size();
+        for (size_type d = 0; d < dimension; ++d)
+        {
+            step(m_offset + d, shape[d] - 1);
+        }
+        if (dimension != 0)
+        {
+            // In column major order the innermost dimension is the first one of the
+            // iterated shape, which is not a view dimension when m_offset != 0.
+            step(l == layout_type::row_major ? m_offset + dimension - 1 : size_type(0));
+        }
     }
 
     template <bool is_const, class CT, class... S>
@@ -2120,175 +2133,66 @@ namespace xt
         return newaxis_count_before<S...>(index + 1) != newaxis_count_before<S...>(index);
     }
 
+    // Moves m_it from the first element of the underlying expression to the
+    // first element of the view.
     template <bool is_const, class CT, class... S>
-    inline void xview_stepper<is_const, CT, S...>::to_end_impl(layout_type l)
+    inline void xview_stepper<is_const, CT, S...>::to_slice_begin()
     {
+        m_index_keeper.fill(0);
         auto func = [](const auto& s) noexcept
         {
-            return xt::value(s, get_size(s) - 1);
-        };
-        auto size_func = [](const auto& s) noexcept
-        {
-            return get_size(s);
+            return xt::value(s, 0);
         };
-
         for (size_type i = 0; i < sizeof...(S); ++i)
         {
             if (!is_newaxis_slice(i))
             {
                 size_type s = apply<size_type>(i, func, p_view->slices());
-                size_type ix = apply<size_type>(i, size_func, p_view->slices());
-                m_index_keeper[i] = ix - size_type(1);
-                size_type index = i - newaxis_count_before<S...>(i);
-                s = p_view->underlying_size(index) - 1 - s;
-                m_it.step_back(index, s);
+                m_it.step(i - newaxis_count_before<S...>(i), s);
             }
         }
-        if (l == layout_type::row_major)
-        {
-            for (size_type i = sizeof...(S); i > 0; --i)
-            {
-                if (!is_newaxis_slice(i - 1))
-                {
-                    m_index_keeper[i - 1]++;
-                    break;
-                }
-            }
-        }
-        else if (l == layout_type::column_major)
-        {
-            for (size_type i = 0; i < sizeof...(S); ++i)
-            {
-                if (!is_newaxis_slice(i))
-                {
-                    m_index_keeper[i]++;
-                    break;
-                }
-            }
-        }
-        else
-        {
-            XTENSOR_THROW(std::runtime_error, "Iteration only allowed in row or column major.");
-        }
     }
 
     template <bool is_const, class CT, class... S>
-    template <class F>
-    void xview_stepper<is_const, CT, S...>::common_step_forward(size_type dim, F f)
+    inline void xview_stepper<is_const, CT, S...>::common_reset(size_type dim, bool backwards)
     {
-        if (dim >= m_offset)
+        if (dim < m_offset)
         {
-            auto func = [&dim, this](const auto& s) noexcept
-            {
-                return step_size(s, this->m_index_keeper[dim]++, 1);
-            };
-            size_type index = integral_skip<S...>(dim);
-            if (!is_newaxis_slice(index))
-            {
-                size_type step_size = index < sizeof...(S) ? apply<size_type>(index, func, p_view->slices())
-                                                           : 1;
-                index -= newaxis_count_before<S...>(index);
-                f(index, step_size);
-            }
+            return;
         }
-    }
-
-    template <bool is_const, class CT, class... S>
-    template <class F>
-    void xview_stepper<is_const, CT, S...>::common_step_forward(size_type dim, size_type n, F f)
-    {
-        if (dim >= m_offset)
+        const size_type index = integral_skip<S...>(dim - m_offset);
+        if (index >= sizeof...(S))
         {
-            auto func = [&dim, &n, this](const auto& s) noexcept
+            // Not covered by a slice: the view spans this whole dimension of the
+            // underlying expression, whose stepper knows its extent.
+            const size_type underlying_dim = index - newaxis_count<S...>();
+            if (backwards)
             {
-                auto st_size = step_size(s, this->m_index_keeper[dim], n);
-                this->m_index_keeper[dim] += n;
-                return size_type(st_size);
-            };
-
-            size_type index = integral_skip<S...>(dim);
-            if (!is_newaxis_slice(index))
-            {
-                size_type step_size = index < sizeof...(S) ? apply<size_type>(index, func, p_view->slices())
-                                                           : n;
-                index -= newaxis_count_before<S...>(index);
-                f(index, step_size);
+                m_it.reset_back(underlying_dim);
             }
-        }
-    }
-
-    template <bool is_const, class CT, class... S>
-    template <class F>
-    void xview_stepper<is_const, CT, S...>::common_step_backward(size_type dim, F f)
-    {
-        if (dim >= m_offset)
-        {
-            auto func = [&dim, this](const auto& s) noexcept
-            {
-                this->m_index_keeper[dim]--;
-                return step_size(s, this->m_index_keeper[dim], 1);
-            };
-            size_type index = integral_skip<S...>(dim);
-            if (!is_newaxis_slice(index))
+            else
             {
-                size_type step_size = index < sizeof...(S) ? apply<size_type>(index, func, p_view->slices())
-                                                           : 1;
-                index -= newaxis_count_before<S...>(index);
-                f(index, step_size);
+                m_it.reset(underlying_dim);
             }
         }
-    }
-
-    template <bool is_const, class CT, class... S>
-    template <class F>
-    void xview_stepper<is_const, CT, S...>::common_step_backward(size_type dim, size_type n, F f)
-    {
-        if (dim >= m_offset)
+        else if (!is_newaxis_slice(index))
         {
-            auto func = [&dim, &n, this](const auto& s) noexcept
+            auto func = [this, index, backwards](const auto& s) noexcept
             {
-                this->m_index_keeper[dim] -= n;
-                return step_size(s, this->m_index_keeper[dim], n);
+                const size_type size = get_size(s);
+                m_index_keeper[index] = backwards ? size - 1 : 0;
+                return xt::value(s, size - 1) - xt::value(s, 0);
             };
-
-            size_type index = integral_skip<S...>(dim);
-            if (!is_newaxis_slice(index))
+            const size_type underlying_dim = index - newaxis_count_before<S...>(index);
+            const size_type n = apply<size_type>(index, func, p_view->slices());
+            if (backwards)
             {
-                size_type step_size = index < sizeof...(S) ? apply<size_type>(index, func, p_view->slices())
-                                                           : n;
-                index -= newaxis_count_before<S...>(index);
-                f(index, step_size);
+                m_it.step(underlying_dim, n);
             }
-        }
-    }
-
-    template <bool is_const, class CT, class... S>
-    template <class F>
-    void xview_stepper<is_const, CT, S...>::common_reset(size_type dim, F f, bool backwards)
-    {
-        auto size_func = [](const auto& s) noexcept
-        {
-            return get_size(s);
-        };
-        auto end_func = [](const auto& s) noexcept
-        {
-            return xt::value(s, get_size(s) - 1) - xt::value(s, 0);
-        };
-
-        size_type index = integral_skip<S...>(dim);
-        if (!is_newaxis_slice(index))
-        {
-            if (dim < m_index_keeper.size())
+            else
             {
-                size_type size = index < sizeof...(S) ? apply<size_type>(index, size_func, p_view->slices())
-                                                      : p_view->shape()[dim];
-                m_index_keeper[dim] = backwards ? size - 1 : 0;
+                m_it.step_back(underlying_dim, n);
             }
-
-            size_type reset_n = index < sizeof...(S) ? apply<size_type>(index, end_func, p_view->slices())
-                                                     : p_view->shape()[dim] - 1;
-            index -= newaxis_count_before<S...>(index);
-            f(index, reset_n);
         }
     }
 }
diff --git a/test/test_xview.cpp b/test/test_xview.cpp
index 8069bb016..69cca47aa 100644
--- a/test/test_xview.cpp
+++ b/test/test_xview.cpp
@@ -1117,16 +1117,14 @@ namespace xt
             ++riter_expv1;
         }
 
-        // TODO: this test cannot pass for views containing squeezes of
-        // keep slices with a single element. This is NOT trivial to fix,
-        // but this use case is not that common.
-        /*auto rciter_expv1 = exp.template rbegin<layout_type::column_major>();
+        auto rciter_expv1 = exp.template rbegin<layout_type::column_major>();
         for (auto iter = v.template rbegin<layout_type::column_major>();
-             iter != v.template rend<layout_type::column_major>(); ++iter)
+             iter != v.template rend<layout_type::column_major>();
+             ++iter)
         {
             EXPECT_EQ(*iter, *rciter_expv1);
             ++rciter_expv1;
-        }*/
+        }
     }
 
     TEST(xview, random_stepper)
@@ -1221,6 +1219,51 @@ namespace xt
         EXPECT_EQ(res, exp);
     }
 
+    TEST(xview, keep_slice_reverse_iteration)
+    {
+        xtensor<int, 2> a = {{0, 1, 2, 3}, {4, 5, 6, 7}, {8, 9, 10, 11}};
+
+        auto v1 = xt::view(a, 1, keep(0, 2, 3));
+        xtensor<int, 1> exp_v1 = {4, 6, 7};
+        test_view_iter(v1, exp_v1);
+
+        auto v2 = xt::view(a, keep(0, 2), 1);
+        xtensor<int, 1> exp_v2 = {1, 9};
+        test_view_iter(v2, exp_v2);
+
+        auto v3 = xt::view(a, keep(0, 2));
+        xtensor<int, 2> exp_v3 = {{0, 1, 2, 3}, {8, 9, 10, 11}};
+        test_view_iter(v3, exp_v3);
+
+        auto v4 = xt::view(a, keep(0, 2), xt::all(), xt::newaxis());
+        xtensor<int, 3> exp_v4 = {{{0}, {1}, {2}, {3}}, {{8}, {9}, {10}, {11}}};
+        test_view_iter(v4, exp_v4);
+
+        auto v5 = xt::view(a, xt::newaxis(), drop(1), 2);
+        xtensor<int, 2> exp_v5 = {{2, 10}};
+        test_view_iter(v5, exp_v5);
+    }
+
+    TEST(xview, keep_broadcast_higher_rank)
+    {
+        xtensor<int, 2> a = {{0, 1, 2, 3}, {4, 5, 6, 7}, {8, 9, 10, 11}};
+
+        auto v = xt::view(a, keep(0, 2), keep(1, 3));
+        xtensor<int, 3> exp = {{{1, 3}, {9, 11}}, {{1, 3}, {9, 11}}};
+        xtensor<int, 3> z = xt::zeros<int>({2, 2, 2});
+        xtensor<int, 3> res = z + v;
+        EXPECT_EQ(res, exp);
+        xtensor<int, 3, layout_type::column_major> zc = xt::zeros<int>({2, 2, 2});
+        xtensor<int, 3, layout_type::column_major> res_cm = zc + v;
+        EXPECT_EQ(res_cm, exp);
+
+        auto w = xt::view(a, 1, keep(0, 2, 3));
+        xtensor<int, 2> b = {{100, 200, 300}, {400, 500, 600}};
+        xtensor<int, 2> res2 = b + w;
+        xtensor<int, 2> exp2 = {{104, 206, 307}, {404, 506, 607}};
+        EXPECT_EQ(res2, exp2);
+    }
+
     TEST(xview, drop_slice)
     {
         xtensor<double, 3, layout_type::row_major> a = {