diff --git a/docs/src/lib/concrete_binary_operations/cartesian_product.md b/docs/src/lib/concrete_binary_operations/cartesian_product.md
index c7378a2169..f5b7a41359 100644
--- a/docs/src/lib/concrete_binary_operations/cartesian_product.md
+++ b/docs/src/lib/concrete_binary_operations/cartesian_product.md
@@ -13,4 +13,6 @@ CurrentModule = LazySets
 cartesian_product(::VPolytope, ::VPolytope)
 cartesian_product(::LazySet, ::LazySet)
 cartesian_product(::SimpleSparsePolynomialZonotope, ::SimpleSparsePolynomialZonotope)
+cartesian_product(::SparsePolynomialZonotope, ::SparsePolynomialZonotope)
+cartesian_product(::SparsePolynomialZonotope, ::AbstractZonotope)
 ```
diff --git a/src/ConcreteOperations/cartesian_product.jl b/src/ConcreteOperations/cartesian_product.jl
index a56faa547a..e266f4a59b 100644
--- a/src/ConcreteOperations/cartesian_product.jl
+++ b/src/ConcreteOperations/cartesian_product.jl
@@ -216,6 +216,14 @@ Compute the Cartesian product of two simple sparse polynomial zonotopes.
 ### Output
 
 The Cartesian product of `P1` and `P2`.
+
+### Algorithm
+
+This method implements Proposition 3.1.22 in [1].
+
+[1] Kochdumper, Niklas. *Extensions of polynomial zonotopes and their application
+    to verification of cyber-physical systems.* PhD diss., Technische Universität
+    München, 2022.
 """
 function cartesian_product(P1::SimpleSparsePolynomialZonotope,
                            P2::SimpleSparsePolynomialZonotope)
@@ -224,3 +232,64 @@ function cartesian_product(P1::SimpleSparsePolynomialZonotope,
     E = cat(expmat(P1), expmat(P2); dims=(1, 2))
     return SimpleSparsePolynomialZonotope(c, G, E)
 end
+
+"""
+    cartesian_product(P1::SparsePolynomialZonotope, P2::SparsePolynomialZonotope)
+
+Compute the Cartesian product of two sparse polynomial zonotopes.
+
+### Input
+
+- `P1` -- sparse polynomial zonotope
+- `P2` -- sparse polynomial zonotope
+
+### Output
+
+The Cartesian product of `P1` and `P2`.
+
+### Algorithm
+
+This method implements Proposition 3.1.22 in [1].
+
+[1] Kochdumper, Niklas. *Extensions of polynomial zonotopes and their application
+    to verification of cyber-physical systems.* PhD diss., Technische Universität
+    München, 2022.
+"""
+function cartesian_product(P1::SparsePolynomialZonotope, P2::SparsePolynomialZonotope)
+    c = vcat(center(P1), center(P2))
+    G = cat(genmat_dep(P1), genmat_dep(P2); dims=(1, 2))
+    GI = cat(genmat_indep(P1), genmat_indep(P2); dims=(1, 2))
+    E = cat(expmat(P1), expmat(P2); dims=(1, 2))
+    return SparsePolynomialZonotope(c, G, GI, E)
+end
+
+"""
+    cartesian_product(SPZ::SparsePolynomialZonotope, Z::AbstractZonotope)
+
+Compute the Cartesian product of a sparse polynomial zonotope and a zonotopic set.
+
+### Input
+
+- `SPZ` -- sparse polynomial zonotope
+- `Z`   -- zonotopic set
+
+### Output
+
+The Cartesian product of `SPZ` and `Z`.
+
+### Algorithm
+
+This method implements Proposition 3.1.22 in [1].
+
+[1] Kochdumper, Niklas. *Extensions of polynomial zonotopes and their application
+    to verification of cyber-physical systems.* PhD diss., Technische Universität
+    München, 2022.
+"""
+@commutative function cartesian_product(SPZ::SparsePolynomialZonotope, Z::AbstractZonotope)
+    c = vcat(center(SPZ), center(Z))
+    G1 = genmat_dep(SPZ)
+    G = vcat(G1, zeros(eltype(G1), size(G1)))
+    GI = cat(genmat_indep(SPZ), genmat(Z); dims=(1, 2))
+    E = expmat(SPZ)
+    return SparsePolynomialZonotope(c, G, GI, E)
+end
diff --git a/src/ConcreteOperations/exact_sum.jl b/src/ConcreteOperations/exact_sum.jl
index 2625178e30..50146f90f4 100644
--- a/src/ConcreteOperations/exact_sum.jl
+++ b/src/ConcreteOperations/exact_sum.jl
@@ -13,11 +13,20 @@ Compute the exact sum of sparse polyomial zonotopes ``P₁`` and ``P₂``.
 ### Output
 
 A `SparsePolynomialZonotope` representing the exact sum ``P₁ ⊞ P₂``.
+
+### Algorithm
+
+This method implements Proposition 3.1.20 in [1].
+
+[1] Kochdumper, Niklas. *Extensions of polynomial zonotopes and their application
+    to verification of cyber-physical systems.* PhD diss., Technische Universität
+    München, 2022.
 """
 function exact_sum(P1::SparsePolynomialZonotope, P2::SparsePolynomialZonotope)
-    indexvector(P1) == indexvector(P2) || throw(ArgumentError("the exact sum " *
-                                                              "is currently only implemented for sparse polynomial zonotopes with " *
-                                                              "the same index vector"))
+    if indexvector(P1) != indexvector(P2)
+        throw(ArgumentError("the exact sum is currently only implemented for " *
+                            "sparse polynomial zonotopes with the same index vector"))
+    end
 
     c = center(P1) + center(P2)
     G = hcat(genmat_dep(P1), genmat_dep(P2))
diff --git a/test/Sets/SparsePolynomialZonotope.jl b/test/Sets/SparsePolynomialZonotope.jl
index f646f85c35..a807e2e42c 100644
--- a/test/Sets/SparsePolynomialZonotope.jl
+++ b/test/Sets/SparsePolynomialZonotope.jl
@@ -50,6 +50,31 @@ for N in [Float64, Float32, Rational{Int}]
     @test genmat_indep(TPZ) == genmat_indep(TPZ)
     @test expmat(TPZ) == expmat(TPZ)
 
+    CPPZ = cartesian_product(PZ, PZ2)
+    @test center(CPPZ) == N[4, 4, 0, 0]
+    @test genmat_dep(CPPZ) == N[2 1 2 0 0 0;
+                                0 2 2 0 0 0;
+                                0 0 0 2 0 1;
+                                0 0 0 1 2 1]
+    @test genmat_indep(CPPZ) == hcat(N[1, 0, 0, 0])
+    @test expmat(CPPZ) == [1 0 3 0 0 0;
+                           0 1 1 0 0 0;
+                           0 0 0 1 0 1;
+                           0 0 0 0 1 3]
+    Z = overapproximate(PZ2, Zonotope)
+    CPPZ = cartesian_product(PZ, Z)
+    @test center(CPPZ) == N[4, 4, 0, 0]
+    @test genmat_dep(CPPZ) == N[2 1 2;
+                                0 2 2;
+                                0 0 0;
+                                0 0 0]
+    @test genmat_indep(CPPZ) == N[1 0 0 0;
+                                  0 0 0 0;
+                                  0 2 0 1;
+                                  0 1 2 1]
+    @test expmat(CPPZ) == [1 0 3;
+                           0 1 1]
+
     MSPZ = minkowski_sum(PZ, PZ2)
     @test center(MSPZ) == [4, 4]
     @test genmat_dep(MSPZ) == [2 1 2 2 0 1; 0 2 2 1 2 1]