diff --git a/src/PowerManifold.jl b/src/PowerManifold.jl
index 61e45e0a..d343cd60 100644
--- a/src/PowerManifold.jl
+++ b/src/PowerManifold.jl
@@ -225,6 +225,47 @@ function allocation_promotion_function(M::AbstractPowerManifold, f, args::Tuple)
     return allocation_promotion_function(M.manifold, f, args)
 end
 
+"""
+    change_representer(M::AbstractPowerManifold, ::AbstractMetric, p, X)
+
+Since the metric on a power manifold decouples, the change of a representer can be done elementwise
+"""
+change_representer(::AbstractPowerManifold, ::AbstractMetric, ::Any, ::Any)
+
+function change_representer!(M::AbstractPowerManifold, Y, G::AbstractMetric, p, X)
+    rep_size = representation_size(M.manifold)
+    for i in get_iterator(M)
+        change_representer!(
+            M.manifold,
+            _write(M, rep_size, Y, i),
+            G,
+            _read(M, rep_size, p, i),
+            _read(M, rep_size, X, i),
+        )
+    end
+    return Y
+end
+
+"""
+    change_metric(M::AbstractPowerManifold, ::AbstractMetric, p, X)
+
+Since the metric on a power manifold decouples, the change of metric can be done elementwise.
+"""
+change_metric(M::AbstractPowerManifold, ::AbstractMetric, ::Any, ::Any)
+
+function change_metric!(M::AbstractPowerManifold, Y, G::AbstractMetric, p, X)
+    rep_size = representation_size(M.manifold)
+    for i in get_iterator(M)
+        change_metric!(
+            M.manifold,
+            _write(M, rep_size, Y, i),
+            G,
+            _read(M, rep_size, p, i),
+            _read(M, rep_size, X, i),
+        )
+    end
+    return Y
+end
 
 """
     check_point(M::AbstractPowerManifold, p; kwargs...)
@@ -1301,6 +1342,10 @@ function Base.show(
     return nothing
 end
 
+function vector_bundle_transport(::VectorSpaceType, ::PowerManifold)
+    return ParallelTransport()
+end
+
 function vector_transport_direction!(
     M::AbstractPowerManifold,
     Y,
@@ -1456,6 +1501,29 @@ function Base.view(
     return view(p[I...], rep_size_to_colons(rep_size)...)
 end
 
+@doc raw"""
+    Y = Weingarten(M::AbstractPowerManifold, p, X, V)
+    Weingarten!(M::AbstractPowerManifold, Y, p, X, V)
+
+Since the metric decouples, also the computation of the Weingarten map
+``\mathcal W_p`` can be computed elementwise on the single elements of the [`PowerManifold`](@ref) `M`.
+"""
+Weingarten(::AbstractPowerManifold, p, X, V)
+
+function Weingarten!(M::AbstractPowerManifold, Y, p, X, V)
+    rep_size = representation_size(M.manifold)
+    for i in get_iterator(M)
+        Weingarten!(
+            M.manifold,
+            _write(M, rep_size, Y, i),
+            _read(M, rep_size, p, i),
+            _read(M, rep_size, X, i),
+            _read(M, rep_size, V, i),
+        )
+    end
+    return Y
+end
+
 @inline function _write(M::AbstractPowerManifold, rep_size::Tuple, x::AbstractArray, i::Int)
     return _write(M, rep_size, x, (i,))
 end
diff --git a/src/VectorBundle.jl b/src/VectorBundle.jl
index 20281d17..43abc7ac 100644
--- a/src/VectorBundle.jl
+++ b/src/VectorBundle.jl
@@ -344,9 +344,6 @@ function representation_size(B::TangentBundleFibers)
     return representation_size(B.manifold)
 end
 
-function Base.show(io::IO, tpt::TensorProductType)
-    return print(io, "TensorProductType(", join(tpt.spaces, ", "), ")")
-end
 function Base.show(io::IO, vb::VectorBundle)
     return print(io, "VectorBundle($(vb.type.fiber), $(vb.manifold))")
 end
@@ -368,14 +365,6 @@ Determine the vector tranport used for [`exp`](@ref exp(::FiberBundle, ::Any...)
 """
 fiber_bundle_transport(::VectorSpaceType, ::AbstractManifold) = ParallelTransport()
 
-function vector_space_dimension(M::AbstractManifold, V::TensorProductType)
-    dim = 1
-    for space in V.spaces
-        dim *= fiber_dimension(M, space)
-    end
-    return dim
-end
-
 function vector_space_dimension(B::TangentBundleFibers)
     return manifold_dimension(B.manifold)
 end
diff --git a/src/VectorFiber.jl b/src/VectorFiber.jl
index d88c11b1..4e5b1705 100644
--- a/src/VectorFiber.jl
+++ b/src/VectorFiber.jl
@@ -1,16 +1,4 @@
 
-"""
-    TensorProductType(spaces::VectorSpaceType...)
-
-Vector space type corresponding to the tensor product of given vector space
-types.
-"""
-struct TensorProductType{TS<:Tuple} <: VectorSpaceType
-    spaces::TS
-end
-
-TensorProductType(spaces::VectorSpaceType...) = TensorProductType{typeof(spaces)}(spaces)
-
 """
     VectorSpaceFiberType{TVS<:VectorSpaceType}
 
diff --git a/test/power.jl b/test/power.jl
index b5494922..f8b32f4b 100644
--- a/test/power.jl
+++ b/test/power.jl
@@ -5,6 +5,8 @@ using StaticArrays
 using LinearAlgebra
 using Random
 
+include("test_manifolds.jl")
+
 power_array_wrapper(::Type{NestedPowerRepresentation}, ::Int) = identity
 power_array_wrapper(::Type{NestedReplacingPowerRepresentation}, i::Int) = SVector{i}
 
@@ -293,4 +295,23 @@ struct TestArrayRepresentation <: AbstractPowerRepresentation end
                   rand(MersenneTwister(123), N; vector_at = p)
         end
     end
+
+    @testset "metric conversion" begin
+        M = TestSPD(3)
+        N = PowerManifold(M, NestedPowerRepresentation(), 2)
+        e = EuclideanMetric()
+        p = [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1]
+        q = [2.0 0.0 0.0; 0.0 2.0 0.0; 0.0 0.0 1]
+        P = [p, q]
+        X = [log(M, p, q), log(M, q, p)]
+        Y = change_metric(N, e, P, X)
+        Yc = [change_metric(M, e, p, log(M, p, q)), change_metric(M, e, q, log(M, q, p))]
+        @test norm(N, P, Y .- Yc) ≈ 0
+        Z = change_representer(N, e, P, X)
+        Zc = [
+            change_representer(M, e, p, log(M, p, q)),
+            change_representer(M, e, q, log(M, q, p)),
+        ]
+        @test norm(N, P, Z .- Zc) ≈ 0
+    end
 end
diff --git a/test/test_manifolds.jl b/test/test_manifolds.jl
new file mode 100644
index 00000000..571aab9e
--- /dev/null
+++ b/test/test_manifolds.jl
@@ -0,0 +1,226 @@
+using ManifoldsBase: ℝ, ℂ, DefaultManifold, RealNumbers, EuclideanMetric
+using LinearAlgebra
+
+# minimal sphere implementation for testing more complicated manifolds
+struct TestSphere{N,𝔽} <: AbstractManifold{𝔽} end
+TestSphere(N::Int, 𝔽 = ℝ) = TestSphere{N,𝔽}()
+
+ManifoldsBase.representation_size(::TestSphere{N}) where {N} = (N + 1,)
+
+function ManifoldsBase.change_representer!(
+    M::TestSphere,
+    Y,
+    ::ManifoldsBase.EuclideanMetric,
+    p,
+    X,
+)
+    return copyto!(M, Y, p, X)
+end
+
+function ManifoldsBase.change_metric!(
+    M::TestSphere,
+    Y,
+    ::ManifoldsBase.EuclideanMetric,
+    p,
+    X,
+)
+    return copyto!(M, Y, p, X)
+end
+
+function ManifoldsBase.check_point(M::TestSphere, p; kwargs...)
+    if !isapprox(norm(p), 1.0; kwargs...)
+        return DomainError(
+            norm(p),
+            "The point $(p) does not lie on the $(M) since its norm is not 1.",
+        )
+    end
+    return nothing
+end
+
+function ManifoldsBase.check_vector(M::TestSphere, p, X; kwargs...)
+    if !isapprox(abs(real(dot(p, X))), 0.0; kwargs...)
+        return DomainError(
+            abs(dot(p, X)),
+            "The vector $(X) is not a tangent vector to $(p) on $(M), since it is not orthogonal in the embedding.",
+        )
+    end
+    return nothing
+end
+
+function ManifoldsBase.exp!(M::TestSphere, q, p, X)
+    return exp!(M, q, p, X, one(number_eltype(X)))
+end
+function ManifoldsBase.exp!(::TestSphere, q, p, X, t::Number)
+    θ = abs(t) * norm(X)
+    if θ == 0
+        copyto!(q, p)
+    else
+        X_scale = t * sin(θ) / θ
+        q .= p .* cos(θ) .+ X .* X_scale
+    end
+    return q
+end
+
+function ManifoldsBase.get_basis_diagonalizing(
+    M::TestSphere{n},
+    p,
+    B::DiagonalizingOrthonormalBasis{ℝ},
+) where {n}
+    A = zeros(n + 1, n + 1)
+    A[1, :] = transpose(p)
+    A[2, :] = transpose(B.frame_direction)
+    V = nullspace(A)
+    κ = ones(n)
+    if !iszero(B.frame_direction)
+        # if we have a nonzero direction for the geodesic, add it and it gets curvature zero from the tensor
+        V = hcat(B.frame_direction / norm(M, p, B.frame_direction), V)
+        κ[1] = 0 # no curvature along the geodesic direction, if x!=y
+    end
+    T = typeof(similar(B.frame_direction))
+    Ξ = [convert(T, V[:, i]) for i in 1:n]
+    return CachedBasis(B, κ, Ξ)
+end
+
+@inline function nzsign(z, absz = abs(z))
+    psignz = z / absz
+    return ifelse(iszero(absz), one(psignz), psignz)
+end
+
+
+function ManifoldsBase.get_coordinates_orthonormal!(M::TestSphere, Y, p, X, ::RealNumbers)
+    n = manifold_dimension(M)
+    p1 = p[1]
+    cosθ = abs(p1)
+    λ = nzsign(p1, cosθ)
+    pend, Xend = view(p, 2:(n + 1)), view(X, 2:(n + 1))
+    factor = λ * X[1] / (1 + cosθ)
+    Y .= Xend .- pend .* factor
+    return Y
+end
+
+function ManifoldsBase.get_vector_orthonormal!(M::TestSphere, Y, p, X, ::RealNumbers)
+    n = manifold_dimension(M)
+    p1 = p[1]
+    cosθ = abs(p1)
+    λ = nzsign(p1, cosθ)
+    pend = view(p, 2:(n + 1))
+    pX = dot(pend, X)
+    factor = pX / (1 + cosθ)
+    Y[1] = -λ * pX
+    Y[2:(n + 1)] .= X .- pend .* factor
+    return Y
+end
+
+ManifoldsBase.inner(::TestSphere, p, X, Y) = dot(X, Y)
+
+ManifoldsBase.injectivity_radius(::TestSphere) = π
+
+function ManifoldsBase.inverse_retract_project!(M::TestSphere, X, p, q)
+    X .= q .- p
+    project!(M, X, p, X)
+    return X
+end
+
+ManifoldsBase.is_flat(M::TestSphere) = manifold_dimension(M) == 1
+
+function ManifoldsBase.log!(::TestSphere, X, p, q)
+    cosθ = clamp(real(dot(p, q)), -1, 1)
+    X .= q .- p .* cosθ
+    nrmX = norm(X)
+    if nrmX > 0
+        X .*= acos(cosθ) / nrmX
+    end
+    return X
+end
+
+ManifoldsBase.manifold_dimension(::TestSphere{N}) where {N} = N
+
+function ManifoldsBase.parallel_transport_to!(::TestSphere, Y, p, X, q)
+    m = p .+ q
+    mnorm2 = real(dot(m, m))
+    factor = 2 * real(dot(X, q)) / mnorm2
+    Y .= X .- m .* factor
+    return Y
+end
+
+function ManifoldsBase.project!(::TestSphere, q, p)
+    q .= p ./ norm(p)
+    return q
+end
+
+function ManifoldsBase.project!(::TestSphere, Y, p, X)
+    Y .= X .- p .* real(dot(p, X))
+    return Y
+end
+
+function Random.rand!(M::TestSphere, pX; vector_at = nothing, σ = one(eltype(pX)))
+    return rand!(Random.default_rng(), M, pX; vector_at = vector_at, σ = σ)
+end
+function Random.rand!(
+    rng::AbstractRNG,
+    M::TestSphere,
+    pX;
+    vector_at = nothing,
+    σ = one(eltype(pX)),
+)
+    if vector_at === nothing
+        project!(M, pX, randn(rng, eltype(pX), representation_size(M)))
+    else
+        n = σ * randn(rng, eltype(pX), size(pX)) # Gaussian in embedding
+        project!(M, pX, vector_at, n) #project to TpM (keeps Gaussianness)
+    end
+    return pX
+end
+
+function ManifoldsBase.riemann_tensor!(M::TestSphere, Xresult, p, X, Y, Z)
+    innerZX = inner(M, p, Z, X)
+    innerZY = inner(M, p, Z, Y)
+    Xresult .= innerZY .* X .- innerZX .* Y
+    return Xresult
+end
+
+# from Absil, Mahony, Trumpf, 2013 https://sites.uclouvain.be/absil/2013-01/Weingarten_07PA_techrep.pdf
+function ManifoldsBase.Weingarten!(::TestSphere, Y, p, X, V)
+    return Y .= -X * p' * V
+end
+
+# minimal implementation of SPD (for its nontrivial metric conversion)
+
+struct TestSPD <: AbstractManifold{ℝ}
+    n::Int
+end
+
+function ManifoldsBase.change_representer!(::TestSPD, Y, ::EuclideanMetric, p, X)
+    Y .= p * X * p
+    return Y
+end
+
+function ManifoldsBase.change_metric!(::TestSPD, Y, ::EuclideanMetric, p, X)
+    Y .= p * X
+    return Y
+end
+
+function ManifoldsBase.inner(::TestSPD, p, X, Y)
+    F = cholesky(Symmetric(convert(AbstractMatrix, p)))
+    return dot((F \ Symmetric(X)), (Symmetric(Y) / F))
+end
+
+function spd_sqrt_and_sqrt_inv(p::AbstractMatrix)
+    e = eigen(Symmetric(p))
+    U = e.vectors
+    S = max.(e.values, floatmin(eltype(e.values)))
+    Ssqrt = Diagonal(sqrt.(S))
+    SsqrtInv = Diagonal(1 ./ sqrt.(S))
+    return (Symmetric(U * Ssqrt * transpose(U)), Symmetric(U * SsqrtInv * transpose(U)))
+end
+
+function ManifoldsBase.log!(::TestSPD, X, p, q)
+    (p_sqrt, p_sqrt_inv) = spd_sqrt_and_sqrt_inv(p)
+    T = Symmetric(p_sqrt_inv * convert(AbstractMatrix, q) * p_sqrt_inv)
+    e2 = eigen(T)
+    Se = Diagonal(log.(max.(e2.values, eps())))
+    pUe = p_sqrt * e2.vectors
+    return mul!(X, pUe, Se * transpose(pUe))
+end
+
+ManifoldsBase.representation_size(M::TestSPD) = (M.n, M.n)
diff --git a/test/test_sphere.jl b/test/test_sphere.jl
index 7c025347..931a2e00 100644
--- a/test/test_sphere.jl
+++ b/test/test_sphere.jl
@@ -3,188 +3,7 @@ using ManifoldsBase
 using ManifoldsBase: ℝ, ℂ, DefaultManifold, RealNumbers
 using Test
 
-# minimal sphere implementation for testing more complicated manifolds
-struct TestSphere{N,𝔽} <: AbstractManifold{𝔽} end
-TestSphere(N::Int, 𝔽 = ℝ) = TestSphere{N,𝔽}()
-
-ManifoldsBase.representation_size(::TestSphere{N}) where {N} = (N + 1,)
-
-function ManifoldsBase.change_representer!(
-    M::TestSphere,
-    Y,
-    ::ManifoldsBase.EuclideanMetric,
-    p,
-    X,
-)
-    return copyto!(M, Y, p, X)
-end
-
-function ManifoldsBase.change_metric!(
-    M::TestSphere,
-    Y,
-    ::ManifoldsBase.EuclideanMetric,
-    p,
-    X,
-)
-    return copyto!(M, Y, p, X)
-end
-
-function ManifoldsBase.check_point(M::TestSphere, p; kwargs...)
-    if !isapprox(norm(p), 1.0; kwargs...)
-        return DomainError(
-            norm(p),
-            "The point $(p) does not lie on the $(M) since its norm is not 1.",
-        )
-    end
-    return nothing
-end
-
-function ManifoldsBase.check_vector(M::TestSphere, p, X; kwargs...)
-    if !isapprox(abs(real(dot(p, X))), 0.0; kwargs...)
-        return DomainError(
-            abs(dot(p, X)),
-            "The vector $(X) is not a tangent vector to $(p) on $(M), since it is not orthogonal in the embedding.",
-        )
-    end
-    return nothing
-end
-
-function ManifoldsBase.exp!(M::TestSphere, q, p, X)
-    return exp!(M, q, p, X, one(number_eltype(X)))
-end
-function ManifoldsBase.exp!(::TestSphere, q, p, X, t::Number)
-    θ = abs(t) * norm(X)
-    if θ == 0
-        copyto!(q, p)
-    else
-        X_scale = t * sin(θ) / θ
-        q .= p .* cos(θ) .+ X .* X_scale
-    end
-    return q
-end
-
-function ManifoldsBase.get_basis_diagonalizing(
-    M::TestSphere{n},
-    p,
-    B::DiagonalizingOrthonormalBasis{ℝ},
-) where {n}
-    A = zeros(n + 1, n + 1)
-    A[1, :] = transpose(p)
-    A[2, :] = transpose(B.frame_direction)
-    V = nullspace(A)
-    κ = ones(n)
-    if !iszero(B.frame_direction)
-        # if we have a nonzero direction for the geodesic, add it and it gets curvature zero from the tensor
-        V = hcat(B.frame_direction / norm(M, p, B.frame_direction), V)
-        κ[1] = 0 # no curvature along the geodesic direction, if x!=y
-    end
-    T = typeof(similar(B.frame_direction))
-    Ξ = [convert(T, V[:, i]) for i in 1:n]
-    return CachedBasis(B, κ, Ξ)
-end
-
-@inline function nzsign(z, absz = abs(z))
-    psignz = z / absz
-    return ifelse(iszero(absz), one(psignz), psignz)
-end
-
-
-function ManifoldsBase.get_coordinates_orthonormal!(M::TestSphere, Y, p, X, ::RealNumbers)
-    n = manifold_dimension(M)
-    p1 = p[1]
-    cosθ = abs(p1)
-    λ = nzsign(p1, cosθ)
-    pend, Xend = view(p, 2:(n + 1)), view(X, 2:(n + 1))
-    factor = λ * X[1] / (1 + cosθ)
-    Y .= Xend .- pend .* factor
-    return Y
-end
-
-function ManifoldsBase.get_vector_orthonormal!(M::TestSphere, Y, p, X, ::RealNumbers)
-    n = manifold_dimension(M)
-    p1 = p[1]
-    cosθ = abs(p1)
-    λ = nzsign(p1, cosθ)
-    pend = view(p, 2:(n + 1))
-    pX = dot(pend, X)
-    factor = pX / (1 + cosθ)
-    Y[1] = -λ * pX
-    Y[2:(n + 1)] .= X .- pend .* factor
-    return Y
-end
-
-ManifoldsBase.inner(::TestSphere, p, X, Y) = dot(X, Y)
-
-ManifoldsBase.injectivity_radius(::TestSphere) = π
-
-function ManifoldsBase.inverse_retract_project!(M::TestSphere, X, p, q)
-    X .= q .- p
-    project!(M, X, p, X)
-    return X
-end
-
-ManifoldsBase.is_flat(M::TestSphere) = manifold_dimension(M) == 1
-
-function ManifoldsBase.log!(::TestSphere, X, p, q)
-    cosθ = clamp(real(dot(p, q)), -1, 1)
-    X .= q .- p .* cosθ
-    nrmX = norm(X)
-    if nrmX > 0
-        X .*= acos(cosθ) / nrmX
-    end
-    return X
-end
-
-ManifoldsBase.manifold_dimension(::TestSphere{N}) where {N} = N
-
-function ManifoldsBase.parallel_transport_to!(::TestSphere, Y, p, X, q)
-    m = p .+ q
-    mnorm2 = real(dot(m, m))
-    factor = 2 * real(dot(X, q)) / mnorm2
-    Y .= X .- m .* factor
-    return Y
-end
-
-function ManifoldsBase.project!(::TestSphere, q, p)
-    q .= p ./ norm(p)
-    return q
-end
-
-function ManifoldsBase.project!(::TestSphere, Y, p, X)
-    Y .= X .- p .* real(dot(p, X))
-    return Y
-end
-
-function Random.rand!(M::TestSphere, pX; vector_at = nothing, σ = one(eltype(pX)))
-    return rand!(Random.default_rng(), M, pX; vector_at = vector_at, σ = σ)
-end
-function Random.rand!(
-    rng::AbstractRNG,
-    M::TestSphere,
-    pX;
-    vector_at = nothing,
-    σ = one(eltype(pX)),
-)
-    if vector_at === nothing
-        project!(M, pX, randn(rng, eltype(pX), representation_size(M)))
-    else
-        n = σ * randn(rng, eltype(pX), size(pX)) # Gaussian in embedding
-        project!(M, pX, vector_at, n) #project to TpM (keeps Gaussianness)
-    end
-    return pX
-end
-
-function ManifoldsBase.riemann_tensor!(M::TestSphere, Xresult, p, X, Y, Z)
-    innerZX = inner(M, p, Z, X)
-    innerZY = inner(M, p, Z, Y)
-    Xresult .= innerZY .* X .- innerZX .* Y
-    return Xresult
-end
-
-# from Absil, Mahony, Trumpf, 2013 https://sites.uclouvain.be/absil/2013-01/Weingarten_07PA_techrep.pdf
-function ManifoldsBase.Weingarten!(::TestSphere, Y, p, X, V)
-    return Y .= -X * p' * V
-end
+include("test_manifolds.jl")
 
 @testset "TestSphere" begin
     @testset "ShootingInverseRetraction" begin
diff --git a/test/vector_bundle.jl b/test/vector_bundle.jl
index 8c5ec11c..5c33d462 100644
--- a/test/vector_bundle.jl
+++ b/test/vector_bundle.jl
@@ -132,16 +132,6 @@ include("test_sphere.jl")
         end
     end
 
-    @testset "tensor product" begin
-        TT = ManifoldsBase.TensorProductType(TangentSpace, TangentSpace)
-        TTs = "TensorProductType(TangentSpace, TangentSpace)"
-        VBF = VectorBundleFibers(TT, M)
-        @test sprint(show, TT) == TTs
-        @test vector_space_dimension(VBF) == 9
-        @test base_manifold(VBF) == M
-        @test sprint(show, VBF) == "VectorBundleFibers($(TTs), $(M))"
-    end
-
     @testset "Weingarten Map" begin
         p0 = [1.0, 0.0, 0.0]
         M = TangentSpaceAtPoint(M, p0)