diff --git a/Project.toml b/Project.toml
index 45875759..52d7da77 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "LinearMaps"
 uuid = "7a12625a-238d-50fd-b39a-03d52299707e"
-version = "3.10.1"
+version = "3.10.2"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
@@ -10,9 +10,13 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [weakdeps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [extensions]
 LinearMapsChainRulesCoreExt = "ChainRulesCore"
+LinearMapsSparseArraysExt = "SparseArrays"
+LinearMapsStatisticsExt = "Statistics"
 
 [compat]
 ChainRulesCore = "1"
diff --git a/ext/LinearMapsSparseArraysExt.jl b/ext/LinearMapsSparseArraysExt.jl
new file mode 100644
index 00000000..f0cbafb5
--- /dev/null
+++ b/ext/LinearMapsSparseArraysExt.jl
@@ -0,0 +1,117 @@
+module LinearMapsSparseArraysExt
+
+import SparseArrays: sparse, blockdiag, SparseMatrixCSC
+using SparseArrays: AbstractSparseMatrix
+
+using LinearAlgebra, LinearMaps
+import LinearMaps: _issymmetric, _ishermitian
+using LinearMaps: WrappedMap, CompositeMap, LinearCombination, ScaledMap, UniformScalingMap,
+    AdjointMap, TransposeMap, BlockMap, BlockDiagonalMap, KroneckerMap, KroneckerSumMap,
+    VecOrMatMap, AbstractVecOrMatOrQ, MapOrVecOrMat
+using LinearMaps: convert_to_lmaps, _tail, _unsafe_mul!
+
+_issymmetric(A::AbstractSparseMatrix) = issymmetric(A)
+_ishermitian(A::AbstractSparseMatrix) = ishermitian(A)
+
+# blockdiagonal concatenation via extension of blockdiag
+
+"""
+    blockdiag(As::Union{LinearMap,AbstractVecOrMatOrQ}...)::BlockDiagonalMap
+
+Construct a (lazy) representation of the diagonal concatenation of the arguments.
+To avoid fallback to the generic `blockdiag`, there must be a `LinearMap`
+object among the first 8 arguments.
+"""
+blockdiag
+
+for k in 1:8 # is 8 sufficient?
+    Is = ntuple(n->:($(Symbol(:A, n))::AbstractVecOrMatOrQ), Val(k-1))
+    # yields (:A1, :A2, :A3, ..., :A(k-1))
+    L = :($(Symbol(:A, k))::LinearMap)
+    # yields :Ak
+    mapargs = ntuple(n ->:($(Symbol(:A, n))), Val(k-1))
+    # yields (:LinearMap(A1), :LinearMap(A2), ..., :LinearMap(A(k-1)))
+
+    @eval function blockdiag($(Is...), $L, As::MapOrVecOrMat...)
+        return BlockDiagonalMap(convert_to_lmaps($(mapargs...))...,
+                                $(Symbol(:A, k)),
+                                convert_to_lmaps(As...)...)
+    end
+end
+
+# conversion to sparse arrays
+# sparse: create sparse matrix representation of LinearMap
+function sparse(A::LinearMap{T}) where {T}
+    M, N = size(A)
+    rowind = Int[]
+    nzval = T[]
+    colptr = Vector{Int}(undef, N+1)
+    v = fill(zero(T), N)
+    Av = Vector{T}(undef, M)
+
+    @inbounds for i in eachindex(v)
+        v[i] = one(T)
+        _unsafe_mul!(Av, A, v)
+        js = findall(!iszero, Av)
+        colptr[i] = length(nzval) + 1
+        if length(js) > 0
+            append!(rowind, js)
+            append!(nzval, Av[js])
+        end
+        v[i] = zero(T)
+    end
+    colptr[N+1] = length(nzval) + 1
+
+    return SparseMatrixCSC(M, N, colptr, rowind, nzval)
+end
+Base.convert(::Type{SparseMatrixCSC}, A::LinearMap) = sparse(A)
+SparseMatrixCSC(A::LinearMap) = sparse(A)
+
+sparse(A::ScaledMap{<:Any, <:Any, <:VecOrMatMap}) =
+    A.λ * sparse(A.lmap.lmap)
+sparse(A::WrappedMap) = sparse(A.lmap)
+Base.convert(::Type{SparseMatrixCSC}, A::WrappedMap) = convert(SparseMatrixCSC, A.lmap)
+for (T, t) in ((:AdjointMap, adjoint), (:TransposeMap, transpose))
+    @eval sparse(A::$T) = $t(convert(SparseMatrixCSC, A.lmap))
+end
+function sparse(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{VecOrMatMap}}})
+    mats = map(A->getfield(A, :lmap), ΣA.maps)
+    return sum(sparse, mats)
+end
+function sparse(AB::CompositeMap{<:Any, <:Tuple{VecOrMatMap, VecOrMatMap}})
+    B, A = AB.maps
+    return sparse(A.lmap)*sparse(B.lmap)
+end
+function sparse(λA::CompositeMap{<:Any, <:Tuple{VecOrMatMap, UniformScalingMap}})
+    A, J = λA.maps
+    return J.λ*sparse(A.lmap)
+end
+function sparse(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, VecOrMatMap}})
+    J, A = Aλ.maps
+    return sparse(A.lmap)*J.λ
+end
+function sparse(A::BlockMap)
+    return hvcat(
+        A.rows,
+        convert(SparseMatrixCSC, first(A.maps)),
+        convert.(AbstractArray, _tail(A.maps))...
+    )
+end
+function sparse(A::BlockDiagonalMap)
+    return blockdiag(convert.(SparseMatrixCSC, A.maps)...)
+end
+Base.convert(::Type{AbstractMatrix}, A::BlockDiagonalMap) = sparse(A)
+function sparse(A::KroneckerMap)
+    return kron(
+        convert(SparseMatrixCSC, first(A.maps)),
+        convert.(AbstractMatrix, _tail(A.maps))...
+    )
+end
+function sparse(L::KroneckerSumMap)
+    A, B = L.maps
+    IA = sparse(Diagonal(ones(Bool, size(A, 1))))
+    IB = sparse(Diagonal(ones(Bool, size(B, 1))))
+    return kron(convert(AbstractMatrix, A), IB) + kron(IA, convert(AbstractMatrix, B))
+end
+
+end # module LinearMapsSparseArraysExt
diff --git a/ext/LinearMapsStatisticsExt.jl b/ext/LinearMapsStatisticsExt.jl
new file mode 100644
index 00000000..dd4776f6
--- /dev/null
+++ b/ext/LinearMapsStatisticsExt.jl
@@ -0,0 +1,12 @@
+module LinearMapsStatisticsExt
+
+import Statistics: mean
+
+using LinearMaps
+using LinearMaps: LinearMapTupleOrVector, LinearCombination
+
+mean(f::F, maps::LinearMapTupleOrVector) where {F} = sum(f, maps) / length(maps)
+mean(maps::LinearMapTupleOrVector) = mean(identity, maps)
+mean(A::LinearCombination) = mean(A.maps)
+
+end # module ChainRulesCore
diff --git a/src/LinearMaps.jl b/src/LinearMaps.jl
index 417633b1..18fc5c9c 100644
--- a/src/LinearMaps.jl
+++ b/src/LinearMaps.jl
@@ -6,9 +6,6 @@ export ⊗, squarekron, kronsum, ⊕, sumkronsum, khatrirao, facesplitting
 using LinearAlgebra
 using LinearAlgebra: AbstractQ
 import LinearAlgebra: mul!
-using SparseArrays
-
-import Statistics: mean
 
 using Base: require_one_based_indexing
 
@@ -422,6 +419,8 @@ LinearMap{T}(f, args...; kwargs...) where {T} = FunctionMap{T}(f, args...; kwarg
 
 @static if !isdefined(Base, :get_extension)
     include("../ext/LinearMapsChainRulesCoreExt.jl")
+    include("../ext/LinearMapsSparseArraysExt.jl")
+    include("../ext/LinearMapsStatisticsExt.jl")
 end
 
 end # module
diff --git a/src/blockmap.jl b/src/blockmap.jl
index 5771af61..874ee543 100644
--- a/src/blockmap.jl
+++ b/src/blockmap.jl
@@ -505,34 +505,17 @@ for k in 1:8 # is 8 sufficient?
     mapargs = ntuple(n ->:($(Symbol(:A, n))), Val(k-1))
     # yields (:LinearMap(A1), :LinearMap(A2), ..., :LinearMap(A(k-1)))
 
-    @eval begin
-        function SparseArrays.blockdiag($(Is...), $L, As::MapOrVecOrMat...)
+    @eval function Base.cat($(Is...), $L, As::MapOrVecOrMat...; dims::Dims{2})
+        if dims == (1,2)
             return BlockDiagonalMap(convert_to_lmaps($(mapargs...))...,
                                     $(Symbol(:A, k)),
                                     convert_to_lmaps(As...)...)
-        end
-
-        function Base.cat($(Is...), $L, As::MapOrVecOrMat...; dims::Dims{2})
-            if dims == (1,2)
-                return BlockDiagonalMap(convert_to_lmaps($(mapargs...))...,
-                                        $(Symbol(:A, k)),
-                                        convert_to_lmaps(As...)...)
-            else
-                throw(ArgumentError("dims keyword in cat of LinearMaps must be (1,2)"))
-            end
+        else
+            throw(ArgumentError("dims keyword in cat of LinearMaps must be (1,2)"))
         end
     end
 end
 
-"""
-    blockdiag(As::Union{LinearMap,AbstractVecOrMatOrQ}...)::BlockDiagonalMap
-
-Construct a (lazy) representation of the diagonal concatenation of the arguments.
-To avoid fallback to the generic `SparseArrays.blockdiag`, there must be a `LinearMap`
-object among the first 8 arguments.
-"""
-SparseArrays.blockdiag
-
 """
     cat(As::Union{LinearMap,AbstractVecOrMatOrQ}...; dims=(1,2))::BlockDiagonalMap
 
diff --git a/src/conversion.jl b/src/conversion.jl
index c96c41c4..d3677206 100644
--- a/src/conversion.jl
+++ b/src/conversion.jl
@@ -19,39 +19,8 @@ Base.convert(::Type{Array}, A::LinearMap) = convert(Matrix, A)
 Base.convert(::Type{AbstractMatrix}, A::LinearMap) = AbstractMatrix(A)
 Base.convert(::Type{AbstractArray}, A::LinearMap) = convert(AbstractMatrix, A)
 
-# sparse: create sparse matrix representation of LinearMap
-function SparseArrays.sparse(A::LinearMap{T}) where {T}
-    M, N = size(A)
-    rowind = Int[]
-    nzval = T[]
-    colptr = Vector{Int}(undef, N+1)
-    v = fill(zero(T), N)
-    Av = Vector{T}(undef, M)
-
-    @inbounds for i in eachindex(v)
-        v[i] = one(T)
-        _unsafe_mul!(Av, A, v)
-        js = findall(!iszero, Av)
-        colptr[i] = length(nzval) + 1
-        if length(js) > 0
-            append!(rowind, js)
-            append!(nzval, Av[js])
-        end
-        v[i] = zero(T)
-    end
-    colptr[N+1] = length(nzval) + 1
-
-    return SparseMatrixCSC(M, N, colptr, rowind, nzval)
-end
-Base.convert(::Type{SparseMatrixCSC}, A::LinearMap) = sparse(A)
-SparseArrays.SparseMatrixCSC(A::LinearMap) = sparse(A)
-
 # special cases
 
-# ScaledMap
-SparseArrays.sparse(A::ScaledMap{<:Any, <:Any, <:VecOrMatMap}) =
-    A.λ * sparse(A.lmap.lmap)
-
 # UniformScalingMap
 Base.convert(::Type{AbstractMatrix}, J::UniformScalingMap) = Diagonal(fill(J.λ, J.M))
 
@@ -61,19 +30,10 @@ Base.convert(::Type{T}, A::WrappedMap) where {T<:Matrix} = convert(T, A.lmap)
 Base.Matrix{T}(A::VectorMap) where {T} = copyto!(Matrix{eltype(T)}(undef, size(A)), A.lmap)
 Base.convert(::Type{T}, A::VectorMap) where {T<:Matrix} = T(A)
 Base.convert(::Type{AbstractMatrix}, A::WrappedMap) = convert(AbstractMatrix, A.lmap)
-SparseArrays.sparse(A::WrappedMap) = sparse(A.lmap)
-Base.convert(::Type{SparseMatrixCSC}, A::WrappedMap) = convert(SparseMatrixCSC, A.lmap)
 
 # TransposeMap & AdjointMap
 for (T, t) in ((AdjointMap, adjoint), (TransposeMap, transpose))
     @eval Base.convert(::Type{AbstractMatrix}, A::$T) = $t(convert(AbstractMatrix, A.lmap))
-    @eval SparseArrays.sparse(A::$T) = $t(convert(SparseMatrixCSC, A.lmap))
-end
-
-# LinearCombination
-function SparseArrays.sparse(ΣA::LinearCombination{<:Any, <:Tuple{Vararg{VecOrMatMap}}})
-    mats = map(A->getfield(A, :lmap), ΣA.maps)
-    return sum(sparse, mats)
 end
 
 # CompositeMap
@@ -99,50 +59,19 @@ function Base.Matrix{T}(AB::CompositeMap{<:Any, <:Tuple{VecOrMatMap, VecOrMatMap
     B, A = AB.maps
     return mul!(Matrix{T}(undef, size(AB)), A.lmap, B.lmap)
 end
-function SparseArrays.sparse(AB::CompositeMap{<:Any, <:Tuple{VecOrMatMap, VecOrMatMap}})
-    B, A = AB.maps
-    return sparse(A.lmap)*sparse(B.lmap)
-end
 function Base.Matrix{T}(λA::CompositeMap{<:Any, <:Tuple{VecOrMatMap, UniformScalingMap}}) where {T}
     A, J = λA.maps
     return mul!(Matrix{T}(undef, size(λA)), J.λ, A.lmap)
 end
-function SparseArrays.sparse(λA::CompositeMap{<:Any, <:Tuple{VecOrMatMap, UniformScalingMap}})
-    A, J = λA.maps
-    return J.λ*sparse(A.lmap)
-end
 function Base.Matrix{T}(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, VecOrMatMap}}) where {T}
     J, A = Aλ.maps
     return mul!(Matrix{T}(undef, size(Aλ)), A.lmap, J.λ)
 end
-function SparseArrays.sparse(Aλ::CompositeMap{<:Any, <:Tuple{UniformScalingMap, VecOrMatMap}})
-    J, A = Aλ.maps
-    return sparse(A.lmap)*J.λ
-end
-
-# BlockMap & BlockDiagonalMap
-function SparseArrays.sparse(A::BlockMap)
-    return hvcat(
-        A.rows,
-        convert(SparseMatrixCSC, first(A.maps)),
-        convert.(AbstractArray, _tail(A.maps))...
-    )
-end
-Base.convert(::Type{AbstractMatrix}, A::BlockDiagonalMap) = sparse(A)
-function SparseArrays.sparse(A::BlockDiagonalMap)
-    return blockdiag(convert.(SparseMatrixCSC, A.maps)...)
-end
 
 # KroneckerMap & KroneckerSumMap
 Base.Matrix{T}(A::KroneckerMap) where {T} = kron(convert.(Matrix{T}, A.maps)...)
 Base.convert(::Type{AbstractMatrix}, A::KroneckerMap) =
     kron(convert.(AbstractMatrix, A.maps)...)
-function SparseArrays.sparse(A::KroneckerMap)
-    return kron(
-        convert(SparseMatrixCSC, first(A.maps)),
-        convert.(AbstractMatrix, _tail(A.maps))...
-    )
-end
 
 function Base.Matrix{T}(L::KroneckerSumMap) where {T}
     A, B = L.maps
@@ -156,12 +85,6 @@ function Base.convert(::Type{AbstractMatrix}, L::KroneckerSumMap)
     IB = Diagonal(ones(Bool, size(B, 1)))
     return kron(convert(AbstractMatrix, A), IB) + kron(IA, convert(AbstractMatrix, B))
 end
-function SparseArrays.sparse(L::KroneckerSumMap)
-    A, B = L.maps
-    IA = sparse(Diagonal(ones(Bool, size(A, 1))))
-    IB = sparse(Diagonal(ones(Bool, size(B, 1))))
-    return kron(convert(AbstractMatrix, A), IB) + kron(IA, convert(AbstractMatrix, B))
-end
 
 # FillMap
 Base.Matrix{T}(A::FillMap) where {T} = fill(T(A.λ), size(A))
diff --git a/src/linearcombination.jl b/src/linearcombination.jl
index 4d19b38c..19bceda3 100644
--- a/src/linearcombination.jl
+++ b/src/linearcombination.jl
@@ -21,10 +21,6 @@ Base.mapreduce(::typeof(identity), ::typeof(Base.add_sum), maps::LinearMapTupleO
 Base.mapreduce(::typeof(identity), ::typeof(Base.add_sum), maps::AbstractVector{<:LinearMap{T}}) where {T} =
     LinearCombination{T}(maps)
 
-mean(f::F, maps::LinearMapTupleOrVector) where {F} = sum(f, maps) / length(maps)
-mean(maps::LinearMapTupleOrVector) = mean(identity, maps)
-mean(A::LinearCombination) = mean(A.maps)
-
 MulStyle(A::LinearCombination) = MulStyle(A.maps...)
 
 # basic methods
diff --git a/src/wrappedmap.jl b/src/wrappedmap.jl
index bc14056a..abe9bafa 100644
--- a/src/wrappedmap.jl
+++ b/src/wrappedmap.jl
@@ -24,14 +24,12 @@ WrappedMap(lmap::MapOrVecOrMat{T}; kwargs...) where {T} = WrappedMap{T}(lmap; kw
 # cheap property checks (usually by type)
 _issymmetric(A::AbstractMatrix) = false
 _issymmetric(A::AbstractQ) = false
-_issymmetric(A::AbstractSparseMatrix) = issymmetric(A)
 _issymmetric(A::LinearMap) = issymmetric(A)
 _issymmetric(A::LinearAlgebra.RealHermSymComplexSym) = issymmetric(A)
 _issymmetric(A::Union{Bidiagonal,Diagonal,SymTridiagonal,Tridiagonal}) = issymmetric(A)
 
 _ishermitian(A::AbstractMatrix) = false
 _ishermitian(A::AbstractQ) = false
-_ishermitian(A::AbstractSparseMatrix) = ishermitian(A)
 _ishermitian(A::LinearMap) = ishermitian(A)
 _ishermitian(A::LinearAlgebra.RealHermSymComplexHerm) = ishermitian(A)
 _ishermitian(A::Union{Bidiagonal,Diagonal,SymTridiagonal,Tridiagonal}) = ishermitian(A)