diff --git a/src/mqt/yaqs/core/data_structures/networks.py b/src/mqt/yaqs/core/data_structures/networks.py
index 6ebf5a4a..0e06399d 100644
--- a/src/mqt/yaqs/core/data_structures/networks.py
+++ b/src/mqt/yaqs/core/data_structures/networks.py
@@ -1738,7 +1738,10 @@ def to_mps(self) -> MPS:
             MPS: An MPS object containing the reshaped tensors.
         """
         converted_tensors: list[NDArray[np.complex128]] = [
-            np.reshape(tensor, (tensor.shape[0] * tensor.shape[1], tensor.shape[2], tensor.shape[3]))
+            np.reshape(
+                tensor,
+                (tensor.shape[0] * tensor.shape[1], tensor.shape[2], tensor.shape[3]),
+            )
             for tensor in self.tensors
         ]
 
@@ -1760,12 +1763,133 @@ def to_matrix(self) -> NDArray[np.complex128]:
         for tensor in self.tensors[1:]:
             mat = oe.contract("abcd, efdg->aebfcg", mat, tensor)
             mat = np.reshape(
-                mat, (mat.shape[0] * mat.shape[1], mat.shape[2] * mat.shape[3], mat.shape[4], mat.shape[5])
+                mat,
+                (
+                    mat.shape[0] * mat.shape[1],
+                    mat.shape[2] * mat.shape[3],
+                    mat.shape[4],
+                    mat.shape[5],
+                ),
             )
 
         # Final left and right bonds should be 1
         return np.squeeze(mat, axis=(2, 3))
 
+    @classmethod
+    def from_matrix(
+        cls,
+        mat: np.ndarray,
+        d: int,
+        max_bond: int | None = None,
+        cutoff: float = 1e-12,
+    ) -> MPO:
+        """Factorize a dense matrix into an MPO with uniform local dimension ``d``.
+
+        Each site has local shape ``(d, d)``.
+        The number of sites ``n`` is inferred from the relation:
+
+            mat.shape = (d**n, d**n)
+
+        Args:
+            mat (np.ndarray):
+                Square matrix of shape ``(d**n, d**n)``.
+            d (int):
+                Physical dimension per site. Must satisfy ``d > 0``.
+            max_bond (int | None):
+                Maximum allowed bond dimension (before truncation).
+            cutoff (float):
+                Singular values ``<= cutoff`` are discarded. By default cutoff=1e-12: all numerically non-zero
+                singular values are included.
+
+        Returns:
+            MPO:
+                An MPO with ``n`` sites, uniform physical dimension ``d`` per site,
+                and bond dimensions determined by SVD truncation.
+
+        Raises:
+            ValueError:
+                If ``d <= 0``;
+                If ``d == 1`` but the matrix is not ``1 x 1``;
+                If the matrix is not square;
+                If ``rows`` is not a power of ``d``;
+                If the inferred number of sites ``n < 1``.
+        """
+        if d <= 0:
+            msg = f"Physical dimension d must be > 0, got d={d}."
+            raise ValueError(msg)
+
+        rows, cols = mat.shape
+
+        if rows != cols:
+            msg = "Matrix must be square for uniform MPO factorization."
+            raise ValueError(msg)
+
+        if d == 1:
+            if rows != 1:
+                msg = "For d == 1 the matrix must be 1x1 since 1**n = 1 for any n."
+                raise ValueError(msg)
+            n = 1
+        else:
+            n_float = np.log(rows) / np.log(d)
+            n = round(n_float)
+
+            if n < 1:
+                msg = f"Inferred chain length n={n} is invalid; matrix dimension {rows} too small for base d={d}."
+                raise ValueError(msg)
+
+            if not np.isclose(n_float, n):
+                msg = f"Matrix dimension {rows} is not a power of d={d}."
+                raise ValueError(msg)
+
+        mat = np.asarray(mat, dtype=np.complex128)
+
+        left_rank = 1
+        rem = mat.reshape(1, rows, cols)
+
+        tensors: list[np.ndarray] = []
+
+        def _truncate(s: np.ndarray) -> int:
+            r = s.size
+            if cutoff > 0.0:
+                r = max(int(np.sum(s > cutoff)), 1)
+            if max_bond is not None:
+                r = min(r, max_bond)
+            return r
+
+        for k in range(n - 1):
+            rest = d ** (n - k - 1)
+
+            rem = rem.reshape(left_rank, d, rest, d, rest)
+            rem_perm = np.transpose(rem, (1, 3, 0, 2, 4))
+            x = rem_perm.reshape(d * d * left_rank, rest * rest)
+
+            u, s, vh = np.linalg.svd(x, full_matrices=False)
+
+            r_keep = _truncate(s)
+
+            u = u[:, :r_keep]
+            s = s[:r_keep]
+            vh = vh[:r_keep, :]
+
+            t_k = u.reshape(d, d, left_rank, r_keep)
+            tensors.append(t_k)
+
+            rem = (s[:, None] * vh).reshape(r_keep, rest, rest)
+            left_rank = r_keep
+
+        rem = rem.reshape(left_rank, d, d)
+        t_last = np.transpose(rem, (1, 2, 0)).reshape(d, d, left_rank, 1)
+        tensors.append(t_last)
+
+        mpo = cls()
+        mpo.tensors = tensors
+        mpo.length = n
+        mpo.physical_dimension = d
+
+        assert mpo.check_if_valid_mpo(), "MPO initialized wrong"
+
+        return mpo
+
     def check_if_valid_mpo(self) -> bool:
         """MPO validity check.
 
diff --git a/tests/core/data_structures/test_networks.py b/tests/core/data_structures/test_networks.py
index 5f7bb904..c979856a 100644
--- a/tests/core/data_structures/test_networks.py
+++ b/tests/core/data_structures/test_networks.py
@@ -204,7 +204,9 @@ def untranspose_block(mpo_tensor: NDArray[np.complex128]) -> NDArray[np.complex1
 
 
 def crandn(
-    size: int | tuple[int, ...], *args: int, seed: np.random.Generator | int | None = None
+    size: int | tuple[int, ...],
+    *args: int,
+    seed: np.random.Generator | int | None = None,
 ) -> NDArray[np.complex128]:
     """Draw random samples from the standard complex normal distribution.
 
@@ -222,7 +224,10 @@ def crandn(
         size = (size,)
     rng = np.random.default_rng(seed)
     # 1 / sqrt(2) is a normalization factor
-    return np.asarray((rng.standard_normal(size) + 1j * rng.standard_normal(size)) / np.sqrt(2), dtype=np.complex128)
+    return np.asarray(
+        (rng.standard_normal(size) + 1j * rng.standard_normal(size)) / np.sqrt(2),
+        dtype=np.complex128,
+    )
 
 
 def random_mps(shapes: list[tuple[int, int, int]], *, normalize: bool = True) -> MPS:
@@ -381,6 +386,59 @@ def test_custom() -> None:
         assert np.allclose(original, created)
 
 
+def test_from_matrix() -> None:
+    """Test that from_matrix() constructs a correct MPO.
+
+    This test constructs a dense Bose-Hubbard Hamiltonian and creates an MPO via from_matrix().
+    It checks:
+    - reconstruction correctness for Bose-Hubbard
+    - random matrices at very large bond dimension
+    - random matrices at moderately truncated bond dimension
+    - all validation error branches (Codecov)
+    """
+    rng = np.random.default_rng()
+
+    length = 5
+    d = 3  # local dimension
+    H = _bose_hubbard_dense(length, d, 0.9, 0.6, 0.2)
+
+    Hmpo = MPO.from_matrix(H, d, 4)
+    assert np.allclose(H, Hmpo.to_matrix())
+
+    H = rng.random((d**length, d**length)) + 1j * rng.random((d**length, d**length))
+    Hmpo = MPO.from_matrix(H, d, 1_000_000)
+    assert np.allclose(H, Hmpo.to_matrix())
+
+    length = 6
+    H = rng.random((d**length, d**length)) + 1j * rng.random((d**length, d**length))
+    Hmpo = MPO.from_matrix(H, d, 728)
+    assert np.max(np.abs(H - Hmpo.to_matrix())) < 1e-2
+
+    mat = np.eye(1)
+    with pytest.raises(ValueError, match="Physical dimension d must be > 0"):
+        MPO.from_matrix(mat, d=0)
+
+    # non-square matrix
+    mat = np.zeros((4, 2))
+    with pytest.raises(ValueError, match="Matrix must be square"):
+        MPO.from_matrix(mat, d=2)
+
+    # d == 1 but matrix not 1x1
+    mat = np.eye(4)
+    with pytest.raises(ValueError, match="1x1"):
+        MPO.from_matrix(mat, d=1)
+
+    # matrix dimension not a power of d
+    mat = np.eye(6)
+    with pytest.raises(ValueError, match="not a power"):
+        MPO.from_matrix(mat, d=2)
+
+    # inferred n < 1 (log(1)/log(100) = 0)
+    mat = np.eye(1)
+    with pytest.raises(ValueError, match="invalid"):
+        MPO.from_matrix(mat, d=100)
+
+
 def test_to_mps() -> None:
     """Test converting an MPO to an MPS.
 
@@ -431,7 +489,12 @@ def test_rotate() -> None:
 
     mpo.rotate(conjugate=False)
     for orig, rotated in zip(original_tensors, mpo.tensors, strict=False):
-        assert rotated.shape == (orig.shape[1], orig.shape[0], orig.shape[2], orig.shape[3])
+        assert rotated.shape == (
+            orig.shape[1],
+            orig.shape[0],
+            orig.shape[2],
+            orig.shape[3],
+        )
         np.testing.assert_allclose(rotated, np.transpose(orig, (1, 0, 2, 3)))
 
     mpo.rotate(conjugate=True)
@@ -474,7 +537,12 @@ def test_mps_initialization(state: str) -> None:
     basis_string = "1001"
 
     if state == "basis":
-        mps = MPS(length=length, physical_dimensions=[pdim] * length, state=state, basis_string=basis_string)
+        mps = MPS(
+            length=length,
+            physical_dimensions=[pdim] * length,
+            state=state,
+            basis_string=basis_string,
+        )
     else:
         mps = MPS(length=length, physical_dimensions=[pdim] * length, state=state)
 
@@ -552,12 +620,20 @@ def test_flip_network() -> None:
     t2 = rng.random(size=(pdim, 2, 2)).astype(np.complex128)
     t3 = rng.random(size=(pdim, 2, 1)).astype(np.complex128)
     original_tensors = [t1, t2, t3]
-    mps = MPS(length, tensors=copy.deepcopy(original_tensors), physical_dimensions=[pdim] * length)
+    mps = MPS(
+        length,
+        tensors=copy.deepcopy(original_tensors),
+        physical_dimensions=[pdim] * length,
+    )
 
     mps.flip_network()
     flipped_tensors = mps.tensors
     assert len(flipped_tensors) == length
-    assert flipped_tensors[0].shape == (pdim, original_tensors[2].shape[2], original_tensors[2].shape[1])
+    assert flipped_tensors[0].shape == (
+        pdim,
+        original_tensors[2].shape[2],
+        original_tensors[2].shape[1],
+    )
     mps.flip_network()
     for orig, now in zip(original_tensors, mps.tensors, strict=False):
         assert np.allclose(orig, now)
@@ -870,7 +946,9 @@ def test_convert_to_vector_fidelity() -> None:
 
     # Define the simulation parameters
     sim_params = StrongSimParams(
-        observables=[Observable(Z(), site) for site in range(num_qubits)], get_state=True, show_progress=False
+        observables=[Observable(Z(), site) for site in range(num_qubits)],
+        get_state=True,
+        show_progress=False,
     )
     simulator.run(state, circ, sim_params)
     assert sim_params.output_state is not None
@@ -894,7 +972,9 @@ def test_convert_to_vector_fidelity_long_range() -> None:
 
     # Define the simulation parameters
     sim_params = StrongSimParams(
-        observables=[Observable(Z(), site) for site in range(num_qubits)], get_state=True, show_progress=False
+        observables=[Observable(Z(), site) for site in range(num_qubits)],
+        get_state=True,
+        show_progress=False,
     )
     simulator.run(state, circ, sim_params)
     assert sim_params.output_state is not None
@@ -1243,7 +1323,10 @@ def test_evaluate_observables_meta_validation_errors() -> None:
 
     # Wrong length (entropy expects exactly two adjacent indices)
     sim_bad_len = AnalogSimParams(
-        [Observable(GateLibrary.entropy(), [1])], elapsed_time=0.1, dt=0.1, show_progress=False
+        [Observable(GateLibrary.entropy(), [1])],
+        elapsed_time=0.1,
+        dt=0.1,
+        show_progress=False,
     )
     results_len = np.empty((1, 1), dtype=np.float64)
     with pytest.raises(AssertionError):
@@ -1251,7 +1334,10 @@ def test_evaluate_observables_meta_validation_errors() -> None:
 
     # Non-adjacent Schmidt cut
     sim_non_adj = AnalogSimParams(
-        [Observable(GateLibrary.schmidt_spectrum(), [0, 2])], elapsed_time=0.1, dt=0.1, show_progress=False
+        [Observable(GateLibrary.schmidt_spectrum(), [0, 2])],
+        elapsed_time=0.1,
+        dt=0.1,
+        show_progress=False,
     )
     results_adj = np.empty((1, 1), dtype=object)
     with pytest.raises(AssertionError):
@@ -1322,7 +1408,9 @@ def test_from_pauli_sum_raises_on_nonpositive_length() -> None:
         mpo.from_pauli_sum(terms=[(1.0, "Z0")], length=-5)
 
 
-def test_from_pauli_sum_raises_on_site_index_out_of_bounds(monkeypatch: pytest.MonkeyPatch) -> None:
+def test_from_pauli_sum_raises_on_site_index_out_of_bounds(
+    monkeypatch: pytest.MonkeyPatch,
+) -> None:
     """Pauli-sum MPO validation: parsed site indices outside [0, L-1] must raise."""
     mpo = MPO()
 
@@ -1333,7 +1421,9 @@ def test_from_pauli_sum_raises_on_site_index_out_of_bounds(monkeypatch: pytest.M
         mpo.from_pauli_sum(terms=[(1.0, "Z0")], length=4)
 
 
-def test_from_pauli_sum_raises_on_invalid_local_op_label(monkeypatch: pytest.MonkeyPatch) -> None:
+def test_from_pauli_sum_raises_on_invalid_local_op_label(
+    monkeypatch: pytest.MonkeyPatch,
+) -> None:
     """Pauli-sum MPO validation: parsed local operator labels must be in _VALID."""
     mpo = MPO()
 
@@ -1367,14 +1457,22 @@ def test_compress_raises_on_invalid_directions() -> None:
     """MPO compress input validation: invalid sweep schedule strings must raise."""
     mpo = MPO()
     mpo.tensors = [np.zeros((2, 2, 1, 1), dtype=complex)]
-    with pytest.raises(ValueError, match=r"directions must be one of \{'lr', 'rl', 'lr_rl', 'rl_lr'\}\."):
+    with pytest.raises(
+        ValueError,
+        match=r"directions must be one of \{'lr', 'rl', 'lr_rl', 'rl_lr'\}\.",
+    ):
         mpo.compress(directions="lr,rl")
 
 
-def test_compress_n_sweeps_zero_returns_without_calling_sweeps(monkeypatch: pytest.MonkeyPatch) -> None:
+def test_compress_n_sweeps_zero_returns_without_calling_sweeps(
+    monkeypatch: pytest.MonkeyPatch,
+) -> None:
     """MPO compress control flow: n_sweeps=0 must return without invoking sweeps."""
     mpo = MPO()
-    mpo.tensors = [np.zeros((2, 2, 1, 1), dtype=complex), np.zeros((2, 2, 1, 1), dtype=complex)]
+    mpo.tensors = [
+        np.zeros((2, 2, 1, 1), dtype=complex),
+        np.zeros((2, 2, 1, 1), dtype=complex),
+    ]
 
     called = False
 
@@ -1393,7 +1491,10 @@ def boom(**_kwargs: object) -> None:
 def test_compress_one_sweep_raises_on_invalid_direction() -> None:
     """MPO _compress_one_sweep input validation: direction must be 'lr' or 'rl'."""
     mpo = MPO()
-    mpo.tensors = [np.zeros((2, 2, 1, 1), dtype=complex), np.zeros((2, 2, 1, 1), dtype=complex)]
+    mpo.tensors = [
+        np.zeros((2, 2, 1, 1), dtype=complex),
+        np.zeros((2, 2, 1, 1), dtype=complex),
+    ]
     with pytest.raises(ValueError, match=r"direction must be 'lr' or 'rl'\."):
         mpo._compress_one_sweep(direction="xx", tol=1e-12, max_bond_dim=None)  # noqa: SLF001