diff --git a/NEWS.md b/NEWS.md
index fc0fd55d24f8f..c9c44d5c7d877 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -81,6 +81,7 @@ Standard library changes
   arithmetic to error if the result may be wrapping. Or use a package such as SaferIntegers.jl when
   constructing the range. ([#40382])
 * TCP socket objects now expose `closewrite` functionality and support half-open mode usage ([#40783]).
+* `extrema` now supports `init` keyword argument ([#36265]).
 * Intersect returns a result with the eltype of the type-promoted eltypes of the two inputs ([#41769]).
 * `Iterators.countfrom` now accepts any type that defines `+`. ([#37747])
 
diff --git a/base/compiler/compiler.jl b/base/compiler/compiler.jl
index c265512afcbf6..6a8ea79c2142a 100644
--- a/base/compiler/compiler.jl
+++ b/base/compiler/compiler.jl
@@ -55,6 +55,24 @@ include("operators.jl")
 include("pointer.jl")
 include("refvalue.jl")
 
+# required for bootstrap
+extrema(itr) = extrema(identity, itr)
+function extrema(f, itr)
+    y = iterate(itr)
+    y === nothing && throw(ArgumentError("collection must be non-empty"))
+    (v, s) = y
+    vmin = vmax = f(v)
+    while true
+        y = iterate(itr, s)
+        y === nothing && break
+        (x, s) = y
+        fx = f(x)
+        vmax = max(fx, vmax)
+        vmin = min(fx, vmin)
+    end
+    return (vmin, vmax)
+end
+
 # checked arithmetic
 const checked_add = +
 const checked_sub = -
diff --git a/base/multidimensional.jl b/base/multidimensional.jl
index 11bfcd418905f..c22d6cd922529 100644
--- a/base/multidimensional.jl
+++ b/base/multidimensional.jl
@@ -1709,9 +1709,13 @@ _unique_dims(A::AbstractArray, dims::Colon) = invoke(unique, Tuple{Any}, A)
 end
 
 """
-    extrema(A::AbstractArray; dims) -> Array{Tuple}
+    extrema([f,] A::AbstractArray; dims, [init]) -> Array{Tuple}
 
-Compute the minimum and maximum elements of an array over the given dimensions.
+Compute the minimum and maximum elements of `A` over dimensions `dims`.
+If `f` is provided, return the minimum and maximum elements after applying `f` to them.
+
+!!! compat "Julia 1.2"
+    The `extrema(f, A)` method requires Julia 1.2 or later.
 
 # Examples
 ```jldoctest
@@ -1734,22 +1738,20 @@ julia> extrema(A, dims = (1,2))
  (9, 15)
 ```
 """
-extrema(A::AbstractArray; dims = :) = _extrema_dims(identity, A, dims)
-
-"""
-    extrema(f, A::AbstractArray; dims) -> Array{Tuple}
-
-Compute the minimum and maximum of `f` applied to each element in the given dimensions
-of `A`.
-
-!!! compat "Julia 1.2"
-    This method requires Julia 1.2 or later.
-"""
-extrema(f, A::AbstractArray; dims=:) = _extrema_dims(f, A, dims)
-
-_extrema_dims(f, A::AbstractArray, ::Colon) = _extrema_itr(f, A)
-
-function _extrema_dims(f, A::AbstractArray, dims)
+extrema(f::F, A::AbstractArray; dims=:, init=_InitialValue()) where {F} =
+    _extrema_dims(f, A, dims, init)
+
+_extrema_dims(f::F, A::AbstractArray, ::Colon, init) where {F} =
+    mapreduce(_DupY(f), _extrema_rf, A; init = init)
+_extrema_dims(f::F, A::AbstractArray, ::Colon, ::_InitialValue) where {F} =
+    mapreduce(_DupY(f), _extrema_rf, A)
+# Note: not passing `init = _InitialValue()` since user-defined
+# `reduce`/`foldl` cannot be aware of `Base._InitialValue` that is an
+# internal implementation detail.
+
+_extrema_dims(f::F, A::AbstractArray, dims, init) where {F} =
+    mapreduce(_DupY(f), _extrema_rf, A; dims = dims, init = init)
+function _extrema_dims(f::F, A::AbstractArray, dims, ::_InitialValue) where {F}
     sz = size(A)
     for d in dims
         sz = setindex(sz, 1, d)
diff --git a/base/operators.jl b/base/operators.jl
index 5bfa7387ec9b7..74f90bb58233f 100644
--- a/base/operators.jl
+++ b/base/operators.jl
@@ -504,58 +504,6 @@ julia> minmax('c','b')
 """
 minmax(x,y) = isless(y, x) ? (y, x) : (x, y)
 
-"""
-    extrema(itr) -> Tuple
-
-Compute both the minimum and maximum element in a single pass, and return them as a 2-tuple.
-
-# Examples
-```jldoctest
-julia> extrema(2:10)
-(2, 10)
-
-julia> extrema([9,pi,4.5])
-(3.141592653589793, 9.0)
-```
-"""
-extrema(itr) = _extrema_itr(identity, itr)
-
-"""
-    extrema(f, itr) -> Tuple
-
-Compute both the minimum and maximum of `f` applied to each element in `itr` and return
-them as a 2-tuple. Only one pass is made over `itr`.
-
-!!! compat "Julia 1.2"
-    This method requires Julia 1.2 or later.
-
-# Examples
-```jldoctest
-julia> extrema(sin, 0:π)
-(0.0, 0.9092974268256817)
-```
-"""
-extrema(f, itr) = _extrema_itr(f, itr)
-
-function _extrema_itr(f, itr)
-    y = iterate(itr)
-    y === nothing && throw(ArgumentError("collection must be non-empty"))
-    (v, s) = y
-    vmin = vmax = f(v)
-    while true
-        y = iterate(itr, s)
-        y === nothing && break
-        (x, s) = y
-        fx = f(x)
-        vmax = max(fx, vmax)
-        vmin = min(fx, vmin)
-    end
-    return (vmin, vmax)
-end
-
-extrema(x::Real) = (x, x)
-extrema(f, x::Real) = (y = f(x); (y, y))
-
 ## definitions providing basic traits of arithmetic operators ##
 
 """
diff --git a/base/reduce.jl b/base/reduce.jl
index 735bd2a0cc9b7..23a56c0c5f45f 100644
--- a/base/reduce.jl
+++ b/base/reduce.jl
@@ -604,7 +604,7 @@ julia> prod(1:5; init = 1.0)
 """
 prod(a; kw...) = mapreduce(identity, mul_prod, a; kw...)
 
-## maximum & minimum
+## maximum, minimum, & extrema
 _fast(::typeof(min),x,y) = min(x,y)
 _fast(::typeof(max),x,y) = max(x,y)
 function _fast(::typeof(max), x::AbstractFloat, y::AbstractFloat)
@@ -785,6 +785,75 @@ Inf
 """
 minimum(a; kw...) = mapreduce(identity, min, a; kw...)
 
+"""
+    extrema(itr; [init]) -> (mn, mx)
+
+Compute both the minimum `mn` and maximum `mx` element in a single pass, and return them
+as a 2-tuple.
+
+The value returned for empty `itr` can be specified by `init`. It must be a 2-tuple whose
+first and second elements are neutral elements for `min` and `max` respectively
+(i.e. which are greater/less than or equal to any other element). As a consequence, when
+`itr` is empty the returned `(mn, mx)` tuple will satisfy `mn ≥ mx`. When `init` is
+specified it may be used even for non-empty `itr`.
+
+!!! compat "Julia 1.8"
+    Keyword argument `init` requires Julia 1.8 or later.
+
+# Examples
+```jldoctest
+julia> extrema(2:10)
+(2, 10)
+
+julia> extrema([9,pi,4.5])
+(3.141592653589793, 9.0)
+
+julia> extrema([]; init = (Inf, -Inf))
+(Inf, -Inf)
+```
+"""
+extrema(itr; kw...) = extrema(identity, itr; kw...)
+
+"""
+    extrema(f, itr; [init]) -> (mn, mx)
+
+Compute both the minimum `mn` and maximum `mx` of `f` applied to each element in `itr` and
+return them as a 2-tuple. Only one pass is made over `itr`.
+
+The value returned for empty `itr` can be specified by `init`. It must be a 2-tuple whose
+first and second elements are neutral elements for `min` and `max` respectively
+(i.e. which are greater/less than or equal to any other element). It is used for non-empty
+collections. Note: it implies that, for empty `itr`, the returned value `(mn, mx)` satisfies
+`mn ≥ mx` even though for non-empty `itr` it  satisfies `mn ≤ mx`.  This is a "paradoxical"
+but yet expected result.
+
+!!! compat "Julia 1.2"
+    This method requires Julia 1.2 or later.
+
+!!! compat "Julia 1.8"
+    Keyword argument `init` requires Julia 1.8 or later.
+
+# Examples
+```jldoctest
+julia> extrema(sin, 0:π)
+(0.0, 0.9092974268256817)
+
+julia> extrema(sin, Real[]; init = (1.0, -1.0))  # good, since -1 ≤ sin(::Real) ≤ 1
+(1.0, -1.0)
+```
+"""
+extrema(f, itr; kw...) = mapreduce(_DupY(f), _extrema_rf, itr; kw...)
+
+# Not using closure since `extrema(type, itr)` is a very likely use-case and it's better
+# to avoid type-instability (#23618).
+struct _DupY{F} <: Function
+    f::F
+end
+_DupY(f::Type{T}) where {T} = _DupY{Type{T}}(f)
+@inline (f::_DupY)(x) = (y = f.f(x); (y, y))
+
+@inline _extrema_rf((min1, max1), (min2, max2)) = (min(min1, min2), max(max1, max2))
+
 ## findmax, findmin, argmax & argmin
 
 """
diff --git a/stdlib/SparseArrays/test/higherorderfns.jl b/stdlib/SparseArrays/test/higherorderfns.jl
index 8da605cf6c0c0..59079682938ee 100644
--- a/stdlib/SparseArrays/test/higherorderfns.jl
+++ b/stdlib/SparseArrays/test/higherorderfns.jl
@@ -709,8 +709,8 @@ end
     @test extrema(f, x) == extrema(f, y)
     @test extrema(spzeros(n, n)) == (0.0, 0.0)
     @test extrema(spzeros(n)) == (0.0, 0.0)
-    @test_throws ArgumentError extrema(spzeros(0, 0))
-    @test_throws ArgumentError extrema(spzeros(0))
+    @test_throws "reducing over an empty" extrema(spzeros(0, 0))
+    @test_throws "reducing over an empty" extrema(spzeros(0))
     @test extrema(sparse(ones(n, n))) == (1.0, 1.0)
     @test extrema(sparse(ones(n))) == (1.0, 1.0)
     @test extrema(A; dims=:) == extrema(B; dims=:)
diff --git a/test/reduce.jl b/test/reduce.jl
index 78ac22c13f366..155ce507a0c0d 100644
--- a/test/reduce.jl
+++ b/test/reduce.jl
@@ -244,23 +244,32 @@ prod2(itr) = invoke(prod, Tuple{Any}, itr)
 
 @test_throws "reducing over an empty" maximum(Int[])
 @test_throws "reducing over an empty" minimum(Int[])
+@test_throws "reducing over an empty" extrema(Int[])
 
 @test maximum(Int[]; init=-1) == -1
 @test minimum(Int[]; init=-1) == -1
+@test extrema(Int[]; init=(1, -1)) == (1, -1)
+
+@test maximum(sin, []; init=-1) == -1
+@test minimum(sin, []; init=1) == 1
+@test extrema(sin, []; init=(1, -1)) == (1, -1)
 
 @test maximum(5) == 5
 @test minimum(5) == 5
 @test extrema(5) == (5, 5)
 @test extrema(abs2, 5) == (25, 25)
+@test Core.Compiler.extrema(abs2, 5) == (25, 25)
 
 let x = [4,3,5,2]
     @test maximum(x) == 5
     @test minimum(x) == 2
     @test extrema(x) == (2, 5)
+    @test Core.Compiler.extrema(x) == (2, 5)
 
     @test maximum(abs2, x) == 25
     @test minimum(abs2, x) == 4
     @test extrema(abs2, x) == (4, 25)
+    @test Core.Compiler.extrema(abs2, x) == (4, 25)
 end
 
 @test maximum([-0.,0.]) === 0.0
diff --git a/test/reducedim.jl b/test/reducedim.jl
index f009a2384ca51..83568a67855b2 100644
--- a/test/reducedim.jl
+++ b/test/reducedim.jl
@@ -93,6 +93,10 @@ A = Array{Int}(undef, 0, 3)
 @test_throws "reducing over an empty collection is not allowed" maximum(A; dims=1)
 @test maximum(A; dims=1, init=-1) == reshape([-1,-1,-1], 1, 3)
 
+@test maximum(zeros(0, 2); dims=1, init=-1) == fill(-1, 1, 2)
+@test minimum(zeros(0, 2); dims=1, init=1) == ones(1, 2)
+@test extrema(zeros(0, 2); dims=1, init=(1, -1)) == fill((1, -1), 1, 2)
+
 # Test reduction along first dimension; this is special-cased for
 # size(A, 1) >= 16
 Breduc = rand(64, 3)