diff --git a/src/sage/matrix/matrix2.pyx b/src/sage/matrix/matrix2.pyx
index 7c5cd644d43..edad4534fe0 100644
--- a/src/sage/matrix/matrix2.pyx
+++ b/src/sage/matrix/matrix2.pyx
@@ -12556,6 +12556,16 @@ cdef class Matrix(Matrix1):
             sage: A.cholesky()
             [   1.0    0.0]
             [-1.0*I    1.0]
+
+        Try the trivial case (:trac:`33107`)::
+
+            sage: all( matrix(R,[]).cholesky() == matrix(R,[])
+            ....:      for R in (RR,CC,RDF,CDF,ZZ,QQ,AA,QQbar) )
+            True
+            sage: all( matrix(R,[]).cholesky().is_immutable()
+            ....:      for R in (RR,CC,RDF,CDF,ZZ,QQ,AA,QQbar) )
+            True
+
         """
         cdef Matrix C # output matrix
         C = self.fetch('cholesky')
@@ -12567,6 +12577,14 @@ cdef class Matrix(Matrix1):
         if not self.is_hermitian():
             raise ValueError("matrix is not Hermitian")
 
+        if n == 0:
+            # The trivial matrix has a trivial cholesky decomposition.
+            # We special-case this after is_hermitian() to ensure that
+            # the matrix is square.
+            C = self.__copy__()
+            C.set_immutable()
+            return C
+
         # Use classical=True to ensure that we don't get a permuted L.
         cdef Matrix L  # block_ldlt() results
         cdef list d    # block_ldlt() results
diff --git a/src/sage/matrix/matrix_double_dense.pyx b/src/sage/matrix/matrix_double_dense.pyx
index 30a3641fd44..e974bfff802 100644
--- a/src/sage/matrix/matrix_double_dense.pyx
+++ b/src/sage/matrix/matrix_double_dense.pyx
@@ -3698,7 +3698,9 @@ cdef class Matrix_double_dense(Matrix_dense):
             raise ValueError(msg.format(self.nrows(), self.ncols()))
         if self._nrows == 0:   # special case
             self.cache(cache_posdef, True)
-            return self.__copy__()
+            L = self.__copy__()
+            L.set_immutable()
+            return L
 
         L = self.fetch(cache_cholesky)
         if L is None: