Repository: vprusso/toqito

Issue:
Feature: Compute the spark of a matrix
Define a function and associated unit tests to compute the *[spark](https://en.wikipedia.org/wiki/Spark_(mathematics))* of a matrix.

Note that the spark is NP-hard to compute in general. Here is a candidate function for computing the spark of a matrix:

```py
def spark(matrix: np.ndarray) -> int:
    """Compute the spark of the given matrix."""

    # The spark of a matrix is equal to 1 if and only if the matrix has a zero column.
    if np.any(np.sum(matrix, axis=0) == 0):
        return 1

    # Otherwise, compute the NP-hard spark computation.
    _, n = matrix.shape
    for k in range(1, n + 1):
        for cols in combinations(range(n), k):
            submatrix = matrix[:, cols]
            if is_linearly_dependent(submatrix):
                return k
    # If no dependent set is found, return n + 1
    return n + 1
```

We would also need unit tests as well. Something like the following, for example:

```py
@pytest.mark.parametrize("matrix, expected_result", [
    # Spark is 3 since (no set of 1 or 2 columns that are linearly dependnet)
    # But there is a set of 3 columns that are linearly dependnet.
    (
        np.array([
            [1, 2, 0, 1],
            [1, 2, 0, 2],
            [1, 2, 0, 3],
            [1, 0, -3, 4],
        ]),
        3
    ),
    # Spark of matrix is 1 if there is a column of all zeros.
    (
        np.array([
            [1, 2, 0, 0],
            [1, 2, 0, 0],
            [1, 2, 0, 0],
            [1, 0, -3, 0],
        ]),
        1
    )
])
def test_spark(matrix, expected_result):
    np.testing.assert_allclose(spark(matrix), expected_result)
```

Presumably, this function would belong in `matrix_props/` and would require a more appropriate and thoroughly written docstring.

Hints from issue comments:
@purva-thakre: I did not know spark of a matrix existed!

The first page on Google search for this quantity also has a link to a character in the matrix movie. lol. 

![image](https://github.com/vprusso/toqito/assets/66048318/ae2ac209-78c5-4150-a2b7-8c821f7a97dd)


Solution patch:
diff --git a/docs/intro_tutorial.rst b/docs/intro_tutorial.rst
index 501627a64..62d825e0a 100644
--- a/docs/intro_tutorial.rst
+++ b/docs/intro_tutorial.rst
@@ -93,7 +93,7 @@ In :code:`toqito`, that can be obtained as
     >>> import numpy as np
     >>> e_0, e_1 = basis(2, 0), basis(2, 1)
     >>> u_0 = 1/np.sqrt(2) * (np.kron(e_0, e_0) + np.kron(e_1, e_1))
-    >>> rho_0 = u_0 * u_0.conj().T
+    >>> rho_0 = u_0 @ u_0.conj().T
     >>> rho_0
     array([[0.5, 0. , 0. , 0.5],
            [0. , 0. , 0. , 0. ],
@@ -194,7 +194,7 @@ Any one of the Bell states serve as an example of a pure state
 
     >>> from toqito.states import bell
     >>> from toqito.state_props import is_pure
-    >>> rho = bell(0) * bell(0).conj().T
+    >>> rho = bell(0) @ bell(0).conj().T
     >>> is_pure(rho)
     True
 
@@ -215,7 +215,7 @@ PPT criterion.
 
     >>> from toqito.states import bell
     >>> from toqito.state_props import is_ppt
-    >>> rho = bell(2) * bell(2).conj().T
+    >>> rho = bell(2) @ bell(2).conj().T
     >>> is_ppt(rho)
     False
 
@@ -237,7 +237,7 @@ For instance, the following bound-entangled tile state is found to be entangled
     >>> from toqito.states import tile
     >>> rho = np.identity(9)
     >>> for i in range(5):
-    ...    rho = rho - tile(i) * tile(i).conj().T
+    ...    rho = rho - tile(i) @ tile(i).conj().T
     >>> rho = rho / 4
     >>> is_separable(rho)
     False
diff --git a/docs/refs.bib b/docs/refs.bib
index 5df631a7c..c926fa9a8 100644
--- a/docs/refs.bib
+++ b/docs/refs.bib
@@ -298,6 +298,13 @@ @article{Eldar_2003_SDPApproach
    year={2003},
    month={Feb}, pages={446-456} }
 
+@book{Elad_2010_Sparse,
+   title={Sparse and redundant representations: from theory to applications in signal and image processing},
+   author={Elad, Michael},
+   year={2010},
+   publisher={Springer Science \& Business Media}
+ }
+
 
 #Last name begins with F
 @article{Fawzi_2015_Positive,
diff --git a/docs/tutorials.extended_nonlocal_games.rst b/docs/tutorials.extended_nonlocal_games.rst
index 78fc4eba1..11698da87 100644
--- a/docs/tutorials.extended_nonlocal_games.rst
+++ b/docs/tutorials.extended_nonlocal_games.rst
@@ -291,13 +291,13 @@ arrays where :code:`prob_mat` corresponds to the probability distribution
     >>> bb84_pred_mat = np.zeros([dim, dim, a_out, b_out, a_in, b_in])
     >>>
     >>> # V(0,0|0,0) = |0><0|
-    >>> bb84_pred_mat[:, :, 0, 0, 0, 0] = e_0 * e_0.conj().T
+    >>> bb84_pred_mat[:, :, 0, 0, 0, 0] = e_0 @ e_0.conj().T
     >>> # V(1,1|0,0) = |1><1|
-    >>> bb84_pred_mat[:, :, 1, 1, 0, 0] = e_1 * e_1.conj().T
+    >>> bb84_pred_mat[:, :, 1, 1, 0, 0] = e_1 @ e_1.conj().T
     >>> # V(0,0|1,1) = |+><+|
-    >>> bb84_pred_mat[:, :, 0, 0, 1, 1] = e_p * e_p.conj().T
+    >>> bb84_pred_mat[:, :, 0, 0, 1, 1] = e_p @ e_p.conj().T
     >>> # V(1,1|1,1) = |-><-|
-    >>> bb84_pred_mat[:, :, 1, 1, 1, 1] = e_m * e_m.conj().T
+    >>> bb84_pred_mat[:, :, 1, 1, 1, 1] = e_m @ e_m.conj().T
     >>>
     >>> # The probability matrix encode \pi(0,0) = \pi(1,1) = 1/2
     >>> bb84_prob_mat = 1/2*np.identity(2)
@@ -666,21 +666,21 @@ Taking the description of :math:`G_{MUB}`, we can encode this as follows.
     >>> num_out = 3
     >>> pred_mat = np.zeros([dim, dim, num_out, num_out, num_in, num_in], dtype=complex)
     >>>
-    >>> pred_mat[:, :, 0, 0, 0, 0] = mubs[0][0] * mubs[0][0].conj().T
-    >>> pred_mat[:, :, 1, 1, 0, 0] = mubs[0][1] * mubs[0][1].conj().T
-    >>> pred_mat[:, :, 2, 2, 0, 0] = mubs[0][2] * mubs[0][2].conj().T
+    >>> pred_mat[:, :, 0, 0, 0, 0] = mubs[0][0] @ mubs[0][0].conj().T
+    >>> pred_mat[:, :, 1, 1, 0, 0] = mubs[0][1] @ mubs[0][1].conj().T
+    >>> pred_mat[:, :, 2, 2, 0, 0] = mubs[0][2] @ mubs[0][2].conj().T
     >>>
-    >>> pred_mat[:, :, 0, 0, 1, 1] = mubs[1][0] * mubs[1][0].conj().T
-    >>> pred_mat[:, :, 1, 1, 1, 1] = mubs[1][1] * mubs[1][1].conj().T
-    >>> pred_mat[:, :, 2, 2, 1, 1] = mubs[1][2] * mubs[1][2].conj().T
+    >>> pred_mat[:, :, 0, 0, 1, 1] = mubs[1][0] @ mubs[1][0].conj().T
+    >>> pred_mat[:, :, 1, 1, 1, 1] = mubs[1][1] @ mubs[1][1].conj().T
+    >>> pred_mat[:, :, 2, 2, 1, 1] = mubs[1][2] @ mubs[1][2].conj().T
     >>>
-    >>> pred_mat[:, :, 0, 0, 2, 2] = mubs[2][0] * mubs[2][0].conj().T
-    >>> pred_mat[:, :, 1, 1, 2, 2] = mubs[2][1] * mubs[2][1].conj().T
-    >>> pred_mat[:, :, 2, 2, 2, 2] = mubs[2][2] * mubs[2][2].conj().T
+    >>> pred_mat[:, :, 0, 0, 2, 2] = mubs[2][0] @ mubs[2][0].conj().T
+    >>> pred_mat[:, :, 1, 1, 2, 2] = mubs[2][1] @ mubs[2][1].conj().T
+    >>> pred_mat[:, :, 2, 2, 2, 2] = mubs[2][2] @ mubs[2][2].conj().T
     >>>
-    >>> pred_mat[:, :, 0, 0, 3, 3] = mubs[3][0] * mubs[3][0].conj().T
-    >>> pred_mat[:, :, 1, 1, 3, 3] = mubs[3][1] * mubs[3][1].conj().T
-    >>> pred_mat[:, :, 2, 2, 3, 3] = mubs[3][2] * mubs[3][2].conj().T
+    >>> pred_mat[:, :, 0, 0, 3, 3] = mubs[3][0] @ mubs[3][0].conj().T
+    >>> pred_mat[:, :, 1, 1, 3, 3] = mubs[3][1] @ mubs[3][1].conj().T
+    >>> pred_mat[:, :, 2, 2, 3, 3] = mubs[3][2] @ mubs[3][2].conj().T
 
 Now that we have encoded :math:`G_{MUB}`, we can calculate the unentangled value.
 
diff --git a/docs/tutorials.state_distinguishability.rst b/docs/tutorials.state_distinguishability.rst
index 3b3bd0157..bd4db1808 100644
--- a/docs/tutorials.state_distinguishability.rst
+++ b/docs/tutorials.state_distinguishability.rst
@@ -123,8 +123,8 @@ Using :code:`toqito`, we can calculate this probability directly as follows:
     >>>
     >>> # Define the corresponding density matrices of |0> and |1> 
     >>> # given as |0><0| and |1><1|, respectively.
-    >>> e_00 = e_0 * e_0.conj().T
-    >>> e_11 = e_1 * e_1.conj().T
+    >>> e_00 = e_0 @ e_0.conj().T
+    >>> e_11 = e_1 @ e_1.conj().T
     >>>
     >>> # Define a list of states and a corresponding list of 
     >>> # probabilities with which those states are selected.
@@ -222,10 +222,10 @@ via PPT measurements in the following manner.
     >>> x_4 = np.kron(psi_3, psi_3)
     >>>
     >>> # YDY density matrices:
-    >>> rho_1 = x_1 * x_1.conj().T
-    >>> rho_2 = x_2 * x_2.conj().T
-    >>> rho_3 = x_3 * x_3.conj().T
-    >>> rho_4 = x_4 * x_4.conj().T
+    >>> rho_1 = x_1 @ x_1.conj().T
+    >>> rho_2 = x_2 @ x_2.conj().T
+    >>> rho_3 = x_3 @ x_3.conj().T
+    >>> rho_4 = x_4 @ x_4.conj().T
     >>>
     >>> states = [rho_1, rho_2, rho_3, rho_4]
     >>> probs = [1 / 4, 1 / 4, 1 / 4, 1 / 4]
diff --git a/poetry.lock b/poetry.lock
index 30f5b1d1a..53ae6a79f 100644
--- a/poetry.lock
+++ b/poetry.lock
@@ -1966,46 +1966,53 @@ pbr = ">=2.0.0"
 
 [[package]]
 name = "symengine"
-version = "0.11.0"
+version = "0.13.0"
 description = "Python library providing wrappers to SymEngine"
 optional = false
-python-versions = ">=3.8,<4"
-files = [
-    {file = "symengine-0.11.0-cp310-cp310-macosx_10_9_x86_64.whl", hash = "sha256:a5ab2281fe5b3cbe8fbe4411c4e2e4ce6788e4d653b5d31b002d05f4af13b7c7"},
-    {file = "symengine-0.11.0-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:1f0c83e91c90a0344fdb9f1d71970f09d4fad5715110c9651d6ebc9ac995f97e"},
-    {file = "symengine-0.11.0-cp310-cp310-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:94c3e143bad732f858f5a8740bda640e24b491aa8cea30b2a68517a552060c84"},
-    {file = "symengine-0.11.0-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:8994e8e9b7273ce48e045d2a97e88b034bb99566bb3d0dd2a2ff79f2dccd4765"},
-    {file = "symengine-0.11.0-cp310-cp310-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:36c96584b5f852990e09904bb890623666cdc0bb83c5758f00939b540fcb2dc9"},
-    {file = "symengine-0.11.0-cp310-cp310-win_amd64.whl", hash = "sha256:fcfaf43781ae6e3e0b6b2f1da700e8ee6a0389426d3300925f4611323eb9b616"},
-    {file = "symengine-0.11.0-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:292efdba67ce3465c99ff1e35338857d982583bbf505ae79f111f8f4fa670367"},
-    {file = "symengine-0.11.0-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:7e935e0f5274d1ab456f3ba75a03dea1ee8a422dbd4d5d4f316653310a177962"},
-    {file = "symengine-0.11.0-cp311-cp311-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:6f588199cb79b1c365005d52a3a8bc0967487dce5db77b864ee73735a5f4d2ec"},
-    {file = "symengine-0.11.0-cp311-cp311-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:bb419847f24944fd29e834a342c6d05577dd66b34942f65131d417d70c2e7926"},
-    {file = "symengine-0.11.0-cp311-cp311-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:10fe4f99d95dc222fd7c37217f80f06f15b9e9b78b951de3c5dc0648cbd1572e"},
-    {file = "symengine-0.11.0-cp311-cp311-win_amd64.whl", hash = "sha256:468b42d2626495808b21cc3b73a162d43c757cc97e4549b1911ec40331d3d491"},
-    {file = "symengine-0.11.0-cp312-cp312-macosx_10_9_x86_64.whl", hash = "sha256:07acd4ab7b147b913fdb62444ced626a9b66b9eefef5d3abfc4f65dccc411403"},
-    {file = "symengine-0.11.0-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:5fb1d62479c18cc00f025b53ba5f72a7bce92eabee9fec44256bb08d3c28d851"},
-    {file = "symengine-0.11.0-cp312-cp312-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:ccd28886f3eea78a1b9150b636ae40acf7d27ac79b4cc9161ec46f592e366a47"},
-    {file = "symengine-0.11.0-cp312-cp312-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:517fec27be3741f89e1eea92b9e2f0cdca70fcd77d3fd5c4f8e4f13e6f04e8e4"},
-    {file = "symengine-0.11.0-cp312-cp312-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:54c5bcca900ec3bc6a20c749957310d52132a815bb3fca614fd66a7881626814"},
-    {file = "symengine-0.11.0-cp312-cp312-win_amd64.whl", hash = "sha256:6123f7e287cebeac9168e60eb6a9e9641f457617dbeee7b117fa7606fde84151"},
-    {file = "symengine-0.11.0-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:116e1f644b3505608b55e745d71b80521b2e72042f421cb70c91dbe7c6710655"},
-    {file = "symengine-0.11.0-cp38-cp38-macosx_11_0_arm64.whl", hash = "sha256:1c7913e264c95f44a52f181e4fff98afd9640651de149d0aee42d73b8b89860b"},
-    {file = "symengine-0.11.0-cp38-cp38-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:ce5e0c14e925fb3600d6955b07314e9981ed389240e437f5dd42f08f70a18af8"},
-    {file = "symengine-0.11.0-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:cc6323c43a703fefb0a50d42430e63a590c18b2955cfa8022cf7c0f6d0f32b86"},
-    {file = "symengine-0.11.0-cp38-cp38-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:86d20163ad6c266e67fb9ad27970c7939bc3d5627487c8f0dc24cd26b9028e3f"},
-    {file = "symengine-0.11.0-cp38-cp38-win_amd64.whl", hash = "sha256:1f79feeeb1af286c29fb791537f243798229020a1ae6340b3d21eece3992783d"},
-    {file = "symengine-0.11.0-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:a78f905944ea9b9ef97e1afb355de4972fe84ed4cceadfc84b01366ce1f50d5b"},
-    {file = "symengine-0.11.0-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:33d2dfea163b77b41fa4c56ade5c6ac4a91693502bc282ba2c645591ef94775d"},
-    {file = "symengine-0.11.0-cp39-cp39-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:2a4412da691328a8e6b0622067807ec61ea97b87dee8d94b1eea2fab2eab8d91"},
-    {file = "symengine-0.11.0-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:c42229ba2fbf70975ebebb0ca5114f4633ba4158246bd7ce14ba9cf4a19828f8"},
-    {file = "symengine-0.11.0-cp39-cp39-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:66def111281fed93da17afd0c787abca48f651078bd3f7987ee4bfb997e26a15"},
-    {file = "symengine-0.11.0-cp39-cp39-win_amd64.whl", hash = "sha256:76b3fceb99d6967d207fbaa93a9cd477949093a63bdb4652cb3c601de85d5229"},
-    {file = "symengine-0.11.0-pp39-pypy39_pp73-macosx_10_9_x86_64.whl", hash = "sha256:7cde52f8c3e8eaa1bf7347a20e50445206cfe058815c61a2a7b86701b1690133"},
-    {file = "symengine-0.11.0-pp39-pypy39_pp73-manylinux_2_12_x86_64.manylinux2010_x86_64.whl", hash = "sha256:83da94f8487b54e391367d66cfdd268636811c4c660daa7b7e90ef9d7b659467"},
-    {file = "symengine-0.11.0-pp39-pypy39_pp73-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:11ae18741ee870ad897f4b1e45ad37587ab78c90a195083d87eb3db2ba66b12c"},
-    {file = "symengine-0.11.0-pp39-pypy39_pp73-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:cb8cc1314312af70b8a148611680f4cd91f8ef971c4fd4befa968cf7d6a2bd19"},
-    {file = "symengine-0.11.0.tar.gz", hash = "sha256:0dd30d29b804ebb7251bddec29c38c3b1fc15ea6953a2c57ee758d5f6fcba458"},
+python-versions = "<4,>=3.8"
+files = [
+    {file = "symengine-0.13.0-cp310-cp310-macosx_10_13_x86_64.whl", hash = "sha256:259fd4111c7a70c72bdff5686de1949e8132baeb612eacdaf8837720c6fe449b"},
+    {file = "symengine-0.13.0-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:44f2eb28a1e36db0bbd6679435412f79da9743bf9c1cb3eff25e0c343b7ddd48"},
+    {file = "symengine-0.13.0-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:d141712fa14d9138bd19e64b10392f850c68d88cd7db29f1bda33e32d1095559"},
+    {file = "symengine-0.13.0-cp310-cp310-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:830226d933bfcdb93546e4062541627d9a3bc7a178a63fb16c002eb5c5221938"},
+    {file = "symengine-0.13.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:a08090163819a0bbfa97d64bd2d8dac2c5268147ed9c242799d7f7e8728a6f4e"},
+    {file = "symengine-0.13.0-cp310-cp310-win_amd64.whl", hash = "sha256:1e435dcd8ed25e4c7c21ab1c0376be910efc7f35da76d532367df27b359f0358"},
+    {file = "symengine-0.13.0-cp311-cp311-macosx_10_13_x86_64.whl", hash = "sha256:da0eba7e106095cdce88eb275c8a9d7c4586ad88f229394c53e1184155c00745"},
+    {file = "symengine-0.13.0-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:b0c175f4f895a73a925508af03faf7efd6cad8593256bbdb5346bd996d3ec5c8"},
+    {file = "symengine-0.13.0-cp311-cp311-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:e58d1e2abd08381aa0cf24c88c0e8b7f592df92619b51e32d36835fbd2dd6ae8"},
+    {file = "symengine-0.13.0-cp311-cp311-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:1db745f2c7a3c5e83510cf4decb43201f43552dfb05ad8af9787c89708be9ede"},
+    {file = "symengine-0.13.0-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:b2572c98b09ac284db6ecff63f6170461194dc94c4209afd34c092ec67873d85"},
+    {file = "symengine-0.13.0-cp311-cp311-win_amd64.whl", hash = "sha256:12727f02a2919f005aee48e68e0cbb70cf857b19385857b4d985d1c9b075f620"},
+    {file = "symengine-0.13.0-cp312-cp312-macosx_10_13_x86_64.whl", hash = "sha256:cf91d24f1bfd6d53228593c7804dd106b71b19674d5afc4fa322d516e1793bdd"},
+    {file = "symengine-0.13.0-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:c5615b7eb68890917abd390ebb10434a949165f6064741c1a8cc345fee14e855"},
+    {file = "symengine-0.13.0-cp312-cp312-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:fb92bdf0890de264abaeacbfbdbd4dd7444b94057bd47958d913b662e549ad8a"},
+    {file = "symengine-0.13.0-cp312-cp312-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:7b3bce486fbc0b87970ed1b10ca9d5cafb1fd6b66382fe631261d83592851d7e"},
+    {file = "symengine-0.13.0-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:f7e6bae9cfcdde2775d92fbb0abe3ef04e32f65ebc4c2d164ca33f4da202d4a7"},
+    {file = "symengine-0.13.0-cp312-cp312-win_amd64.whl", hash = "sha256:bced0a1dbdb94737c299384c85ddbad6944ce8dadc334f7bb8dbbd8f6c965807"},
+    {file = "symengine-0.13.0-cp313-cp313-macosx_10_13_x86_64.whl", hash = "sha256:5d34df77971538e4c29f2d8e5ef7f459c2179465e6cdb7dfd48b79b87ecd8f4d"},
+    {file = "symengine-0.13.0-cp313-cp313-macosx_11_0_arm64.whl", hash = "sha256:ab2661d9b18867e7c6edbfa7a74b8b0a2a694bd24aa08003dc3214f77cb9d6f2"},
+    {file = "symengine-0.13.0-cp313-cp313-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:53f27b9013878ee4419d8e853664d8ae4b68419e3f4b9b5b7f503d32bf904755"},
+    {file = "symengine-0.13.0-cp313-cp313-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:27987f75ce08c64f453e2b9b74fec6ffc5ca418c4deca0b75580979d4a4e242a"},
+    {file = "symengine-0.13.0-cp313-cp313-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:e9ea9410330ea15ed4137d7a0a3c43caccacb71490e18036ce5182d08c93baf8"},
+    {file = "symengine-0.13.0-cp313-cp313-win_amd64.whl", hash = "sha256:5031eb7a5c6675d5195bb57f93cc7d9ac5a7a9a826d4ad6f6b2927746ed7e6e6"},
+    {file = "symengine-0.13.0-cp313-cp313t-macosx_10_13_x86_64.whl", hash = "sha256:ce0e5dfb19943bcf3e44a4485bcac4c5533ba3705c63083494eed0b3bf246076"},
+    {file = "symengine-0.13.0-cp313-cp313t-macosx_11_0_arm64.whl", hash = "sha256:c3b77dc54bf1181f6bd3b3338c4e6e5973a8b0fa20a189d15563ef5626e57b04"},
+    {file = "symengine-0.13.0-cp313-cp313t-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:ca7c3f6c168f6f5b06b421833c3d3baae56067a94b671bdffbe09b8e4fefd9be"},
+    {file = "symengine-0.13.0-cp313-cp313t-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:847523de682416811bacb3ad11507e663b3522fbb35cd27184757e9956d0eaf0"},
+    {file = "symengine-0.13.0-cp313-cp313t-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:ec2fc1b7d96426463f0c9011e9fb88459d906477c1baa8a996dde6fb2bfa99d4"},
+    {file = "symengine-0.13.0-cp38-cp38-macosx_10_13_x86_64.whl", hash = "sha256:6e371bb2da3867085779c1c21bbb1c85f9634c76c8a76c08562ea113e3dfcd85"},
+    {file = "symengine-0.13.0-cp38-cp38-macosx_11_0_arm64.whl", hash = "sha256:cf7c62478b19683d54e4d93faa5b89303beae25db0c503a105a70d266dc99fa9"},
+    {file = "symengine-0.13.0-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:fdb21158cf3e2ba87e441f21ecc7724f108b8db17c0fd1880f9f531602bab1f3"},
+    {file = "symengine-0.13.0-cp38-cp38-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:1848d366b359ff69ef5dac148b30ca04c7339a7d3bcab28419d411e68c0cc011"},
+    {file = "symengine-0.13.0-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:242f817e890a0a50d52ed6b2bfd1aad19636a58db700c7995bbe1ceeaebd9d08"},
+    {file = "symengine-0.13.0-cp38-cp38-win_amd64.whl", hash = "sha256:af79cf2b9645fb55216850185987b2e9347db71e42e87b6402e4bbd41710b316"},
+    {file = "symengine-0.13.0-cp39-cp39-macosx_10_13_x86_64.whl", hash = "sha256:17ae9d1c781a60ac48d07fa30f39a2d237f9da95e9e81f6a24b1c16908e9cad2"},
+    {file = "symengine-0.13.0-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:5264157a95f6d09dd044cee6abcbc176e649c487638b7f32199f387f37ad82a5"},
+    {file = "symengine-0.13.0-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:f01a6d0829b4881d5f831ae7848eb0d82b80d8b46b5689f1bf27069672370f75"},
+    {file = "symengine-0.13.0-cp39-cp39-manylinux_2_17_ppc64le.manylinux2014_ppc64le.whl", hash = "sha256:35ec68b11c2df2be1a236d0c028edb5b331909b16666d7a9fe99a4a5810afec7"},
+    {file = "symengine-0.13.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:b4357fed87083e8719fcffd8bd0e7ddd16172e319343362512f681e472ac5668"},
+    {file = "symengine-0.13.0-cp39-cp39-win_amd64.whl", hash = "sha256:a41be31816e5e51e9063bf26de07faf3751de7a133dbbec149632de702a28e18"},
+    {file = "symengine-0.13.0.tar.gz", hash = "sha256:ab83a08897ebf12579702c2b71ba73d4732fb706cc4291d810aedf39c690c14c"},
 ]
 
 [[package]]
diff --git a/setup.py b/setup.py
index 3cac4876b..89e1624fb 100644
--- a/setup.py
+++ b/setup.py
@@ -11,7 +11,6 @@
     "numpy",
     "picos",
     "scipy",
-    "scikit-image",
     "pytest",
     "pytest-cov",
     "qiskit",
diff --git a/toqito/channel_ops/kraus_to_choi.py b/toqito/channel_ops/kraus_to_choi.py
index 3fd973027..fac763797 100644
--- a/toqito/channel_ops/kraus_to_choi.py
+++ b/toqito/channel_ops/kraus_to_choi.py
@@ -63,7 +63,7 @@ def kraus_to_choi(kraus_ops: list[list[np.ndarray]], sys: int = 2) -> np.ndarray
     dim_op_1, dim_op_2 = dim_in
 
     choi_mat = partial_channel(
-        max_entangled(dim_op_1, False, False) * max_entangled(dim_op_2, False, False).conj().T,
+        max_entangled(dim_op_1, False, False) @ max_entangled(dim_op_2, False, False).conj().T,
         kraus_ops,
         sys,
         np.array([[dim_op_1, dim_op_1], [dim_op_2, dim_op_2]]),
diff --git a/toqito/channel_ops/partial_channel.py b/toqito/channel_ops/partial_channel.py
index 57f6fc90e..215f54c65 100644
--- a/toqito/channel_ops/partial_channel.py
+++ b/toqito/channel_ops/partial_channel.py
@@ -181,7 +181,7 @@ def partial_channel(
         psi_c2 = max_entangled(prod_dim_c2, False, False)
 
         phi_map = permute_systems(
-            np.kron(np.kron(psi_r1 * psi_c1.conj().T, phi_map), psi_r2 * psi_c2.conj().T),
+            np.kron(np.kron(psi_r1 @ psi_c1.conj().T, phi_map), psi_r2 @ psi_c2.conj().T),
             [0, 2, 4, 1, 3, 5],
             dim,
         )
diff --git a/toqito/channels/choi.py b/toqito/channels/choi.py
index 8e2d153ea..c6728ddd1 100644
--- a/toqito/channels/choi.py
+++ b/toqito/channels/choi.py
@@ -103,4 +103,4 @@ def choi(a_var: int = 1, b_var: int = 1, c_var: int = 0) -> np.ndarray:
 
     """
     psi = max_entangled(3, False, False)
-    return np.diag([a_var + 1, c_var, b_var, b_var, a_var + 1, c_var, c_var, b_var, a_var + 1]) - psi * psi.conj().T
+    return np.diag([a_var + 1, c_var, b_var, b_var, a_var + 1, c_var, c_var, b_var, a_var + 1]) - psi @ psi.conj().T
diff --git a/toqito/channels/dephasing.py b/toqito/channels/dephasing.py
index f46bb9bc1..019cf3ab7 100644
--- a/toqito/channels/dephasing.py
+++ b/toqito/channels/dephasing.py
@@ -86,4 +86,4 @@ def dephasing(dim: int, param_p: float = 0) -> np.ndarray:
 
     # Gives a sparse non-normalized state.
     psi = max_entangled(dim=dim, is_sparse=False, is_normalized=False)
-    return (1 - param_p) * np.diag(np.diag(psi * psi.conj().T)) + param_p * (psi * psi.conj().T)
+    return (1 - param_p) * np.diag(np.diag(psi @ psi.conj().T)) + param_p * (psi @ psi.conj().T)
diff --git a/toqito/channels/depolarizing.py b/toqito/channels/depolarizing.py
index af687017e..4e6ce8783 100644
--- a/toqito/channels/depolarizing.py
+++ b/toqito/channels/depolarizing.py
@@ -92,4 +92,4 @@ def depolarizing(dim: int, param_p: float = 0) -> np.ndarray:
 
     # Gives a sparse non-normalized state.
     psi = max_entangled(dim=dim, is_sparse=False, is_normalized=False)
-    return (1 - param_p) * np.identity(dim**2) / dim + param_p * (psi * psi.conj().T)
+    return (1 - param_p) * np.identity(dim**2) / dim + param_p * (psi @ psi.conj().T)
diff --git a/toqito/matrix_props/__init__.py b/toqito/matrix_props/__init__.py
index 2647e0f25..d3b5d52a6 100644
--- a/toqito/matrix_props/__init__.py
+++ b/toqito/matrix_props/__init__.py
@@ -30,4 +30,4 @@
 from toqito.matrix_props.is_positive import is_positive
 from toqito.matrix_props.positive_semidefinite_rank import positive_semidefinite_rank
 from toqito.matrix_props.is_stochastic import is_stochastic
-
+from toqito.matrix_props.spark import spark
diff --git a/toqito/matrix_props/is_density.py b/toqito/matrix_props/is_density.py
index fb8c72820..7224994b5 100644
--- a/toqito/matrix_props/is_density.py
+++ b/toqito/matrix_props/is_density.py
@@ -36,7 +36,7 @@ def is_density(mat: np.ndarray) -> bool:
     >>> from toqito.matrix_props import is_density
     >>> from toqito.states import bell
     >>> import numpy as np
-    >>> rho = bell(0) * bell(0).conj().T
+    >>> rho = bell(0) @ bell(0).conj().T
     >>> is_density(rho)
     np.True_
 
diff --git a/toqito/matrix_props/is_unitary.py b/toqito/matrix_props/is_unitary.py
index ce13ee9b8..4d8901e23 100644
--- a/toqito/matrix_props/is_unitary.py
+++ b/toqito/matrix_props/is_unitary.py
@@ -83,5 +83,5 @@ def is_unitary(mat: np.ndarray, rtol: float = 1e-05, atol: float = 1e-08) -> boo
     u_uc_mat = mat @ mat.conj().T
     id_mat = np.eye(len(mat))
 
-    # If U^* * U = I U * U^*, the matrix "U" is unitary.
+    # If U^* @ U = I U @ U^*, the matrix "U" is unitary.
     return np.allclose(uc_u_mat, id_mat, rtol=rtol, atol=atol) and np.allclose(u_uc_mat, id_mat, rtol=rtol, atol=atol)
diff --git a/toqito/matrix_props/majorizes.py b/toqito/matrix_props/majorizes.py
index f9a050e66..89195fcd2 100644
--- a/toqito/matrix_props/majorizes.py
+++ b/toqito/matrix_props/majorizes.py
@@ -36,7 +36,7 @@ def majorizes(a_var: np.ndarray | list[int], b_var: np.ndarray | list[int]) -> b
     >>> from toqito.channels import partial_trace
     >>>
     >>> v_vec = max_entangled(3)
-    >>> rho = v_vec * v_vec.conj().T
+    >>> rho = v_vec @ v_vec.conj().T
     >>> majorizes(partial_trace(rho, [1]), rho)
     False
 
diff --git a/toqito/matrix_props/spark.py b/toqito/matrix_props/spark.py
new file mode 100644
index 000000000..0da71d92c
--- /dev/null
+++ b/toqito/matrix_props/spark.py
@@ -0,0 +1,61 @@
+"""Compute spark of matrix."""
+
+from itertools import combinations
+
+import numpy as np
+
+
+def spark(mat: np.ndarray) -> int:
+    """Compute the spark of a matrix.
+
+    The spark of a matrix A is the smallest number of columns from A that are linearly
+    dependent :cite:`Elad_2010_Sparse`.
+
+    Examples
+    =========
+    >>> import numpy as np
+    >>> from toqito.matrix_props import spark
+    >>> A = np.array([[1, 0, 1, 2],
+    ...               [0, 1, 1, 3],
+    ...               [1, 1, 2, 5]])
+    >>> spark(A)
+    3
+
+    Notes
+    =====
+    - This function only works for 2D NumPy arrays.
+    - If all columns are linearly independent, the function returns n_cols + 1.
+    - The time complexity of this implementation is O(2^n) in the worst case,
+      where n is the number of columns.
+    - For an m x n matrix A with n >= m:
+      * If spark(A) = m + 1, then rank(A) = m (full rank).
+      * spark(A) = 1 if and only if the matrix has a zero column.
+      * spark(A) <= rank(A) + 1.
+
+    References
+    ==========
+    .. bibliography::
+        :filter: docname in docnames
+
+    :param mat: The input matrix as a 2D NumPy array.
+    :return: The spark of the input matrix :code:`mat`.
+    :raises ValueError: If the input is not a 2D NumPy array.
+
+    """
+    if not isinstance(mat, np.ndarray) or mat.ndim != 2:
+        raise ValueError("Input must be a 2D NumPy array")
+
+    m, n_cols = mat.shape
+
+    # Check for zero columns
+    if np.any(np.all(mat == 0, axis=0)):
+        return 1
+
+    for k in range(1, min(m, n_cols) + 1):
+        for cols in combinations(range(n_cols), k):
+            submatrix = mat[:, cols]
+            if np.linalg.matrix_rank(submatrix) < k:
+                return k
+
+    # If all columns are linearly independent
+    return min(m, n_cols) + 1
diff --git a/toqito/matrix_props/trace_norm.py b/toqito/matrix_props/trace_norm.py
index 5ed35285b..d3472025f 100644
--- a/toqito/matrix_props/trace_norm.py
+++ b/toqito/matrix_props/trace_norm.py
@@ -34,7 +34,7 @@ def trace_norm(rho: np.ndarray) -> float:
 
     >>> from toqito.states import bell
     >>> from toqito.matrix_props import trace_norm
-    >>> rho = bell(0) * bell(0).conj().T
+    >>> rho = bell(0) @ bell(0).conj().T
     >>> trace_norm(rho)
     np.float64(0.9999999999999999)
 
diff --git a/toqito/measurement_ops/measure.py b/toqito/measurement_ops/measure.py
index 9dc56d4c0..f6247f821 100644
--- a/toqito/measurement_ops/measure.py
+++ b/toqito/measurement_ops/measure.py
@@ -40,10 +40,10 @@ def measure(measurement: np.ndarray, state: np.ndarray) -> float:
     >>> e_0, e_1 = basis(2, 0), basis(2, 1)
     >>>
     >>> u = 1/np.sqrt(3) * e_0 + np.sqrt(2/3) * e_1
-    >>> rho = u * u.conj().T
+    >>> rho = u @ u.conj().T
     >>>
-    >>> proj_0 = e_0 * e_0.conj().T
-    >>> proj_1 = e_1 * e_1.conj().T
+    >>> proj_0 = e_0 @ e_0.conj().T
+    >>> proj_1 = e_1 @ e_1.conj().T
 
     Then the probability of obtaining outcome :math:`0` is given by
 
@@ -73,4 +73,4 @@ def measure(measurement: np.ndarray, state: np.ndarray) -> float:
              the variable :code:`measurement` to the variable :code:`state`.
 
     """
-    return float(np.trace(measurement.conj().T * state))
+    return float(np.trace(measurement.conj().T @ state))
diff --git a/toqito/nonlocal_games/extended_nonlocal_game.py b/toqito/nonlocal_games/extended_nonlocal_game.py
index d778b62ac..933edf7e2 100644
--- a/toqito/nonlocal_games/extended_nonlocal_game.py
+++ b/toqito/nonlocal_games/extended_nonlocal_game.py
@@ -112,7 +112,7 @@ def unentangled_value(self) -> float:
 
                 rho = cvxpy.Variable((dim_x, dim_y), hermitian=True)
 
-                objective = cvxpy.Maximize(cvxpy.real(cvxpy.trace(p_win.conj().T @ rho)))
+                objective = cvxpy.Maximize(cvxpy.real(cvxpy.trace(p_win.conj().T * rho)))
 
                 constraints = [cvxpy.trace(rho) == 1, rho >> 0]
                 problem = cvxpy.Problem(objective, constraints)
@@ -194,7 +194,7 @@ def nonsignaling_value(self) -> float:
                 for x_in in range(alice_in):
                     for y_in in range(bob_in):
                         p_win += self.prob_mat[x_in, y_in] * cvxpy.trace(
-                            self.pred_mat[:, :, a_out, b_out, x_in, y_in].conj().T @ k_var[a_out, b_out, x_in, y_in]
+                            self.pred_mat[:, :, a_out, b_out, x_in, y_in].conj().T * k_var[a_out, b_out, x_in, y_in]
                         )
 
         objective = cvxpy.Maximize(cvxpy.real(p_win))
@@ -328,7 +328,7 @@ def __optimize_alice(self, bob_povms) -> tuple[dict, float]:
                                 )
                                 .conj()
                                 .T
-                                @ rho[x_ques, a_ans]
+                                * rho[x_ques, a_ans]
                             )
                         if isinstance(
                             bob_povms[y_ques, b_ans],
@@ -343,7 +343,7 @@ def __optimize_alice(self, bob_povms) -> tuple[dict, float]:
                                 )
                                 .conj()
                                 .T
-                                @ rho[x_ques, a_ans]
+                                * rho[x_ques, a_ans]
                             )
         objective = cvxpy.Maximize(cvxpy.real(win))
         constraints = []
@@ -396,7 +396,7 @@ def __optimize_bob(self, rho) -> tuple[dict, float]:
                                     bob_povms[y_ques, b_ans],
                                 )
                             )
-                            @ rho[x_ques, a_ans].value
+                            * rho[x_ques, a_ans].value
                         )
         objective = cvxpy.Maximize(cvxpy.real(win))
 
@@ -455,7 +455,7 @@ def commuting_measurement_value_upper_bound(self, k: int | str = 1) -> float:
                     for y_in in range(bob_in):
                         p_win += self.prob_mat[x_in, y_in] * cvxpy.trace(
                             self.pred_mat[:, :, a_out, b_out, x_in, y_in].conj().T
-                            @ mat[x_in, y_in][
+                            * mat[x_in, y_in][
                                 a_out * referee_dim : (a_out + 1) * referee_dim,
                                 b_out * referee_dim : (b_out + 1) * referee_dim,
                             ]
diff --git a/toqito/nonlocal_games/nonlocal_game.py b/toqito/nonlocal_games/nonlocal_game.py
index 0109fc03c..7e42f0cfd 100644
--- a/toqito/nonlocal_games/nonlocal_game.py
+++ b/toqito/nonlocal_games/nonlocal_game.py
@@ -372,7 +372,7 @@ def __optimize_alice(self, dim, bob_povms) -> tuple[dict, float]:
                             win += (
                                 self.prob_mat[x_ques, y_ques]
                                 * self.pred_mat[a_ans, b_ans, x_ques, y_ques]
-                                * cvxpy.trace(bob_povms[y_ques, b_ans].conj().T @ alice_povms[x_ques, a_ans])
+                                * cvxpy.trace(bob_povms[y_ques, b_ans].conj().T * alice_povms[x_ques, a_ans])
                             )
                         if isinstance(
                             bob_povms[y_ques, b_ans],
@@ -382,7 +382,7 @@ def __optimize_alice(self, dim, bob_povms) -> tuple[dict, float]:
                             win += (
                                 self.prob_mat[x_ques, y_ques]
                                 * self.pred_mat[a_ans, b_ans, x_ques, y_ques]
-                                * cvxpy.trace(bob_povms[y_ques, b_ans].value.conj().T @ alice_povms[x_ques, a_ans])
+                                * cvxpy.trace(bob_povms[y_ques, b_ans].value.conj().T * alice_povms[x_ques, a_ans])
                             )
 
         objective = cvxpy.Maximize(cvxpy.real(win))
@@ -430,7 +430,7 @@ def __optimize_bob(self, dim, alice_povms) -> tuple[dict, float]:
                         win += (
                             self.prob_mat[x_ques, y_ques]
                             * self.pred_mat[a_ans, b_ans, x_ques, y_ques]
-                            * cvxpy.trace(bob_povms[y_ques, b_ans].H @ alice_povms[x_ques, a_ans].value)
+                            * cvxpy.trace(bob_povms[y_ques, b_ans].H * alice_povms[x_ques, a_ans].value)
                         )
 
         objective = cvxpy.Maximize(cvxpy.real(win))
diff --git a/toqito/nonlocal_games/quantum_hedging.py b/toqito/nonlocal_games/quantum_hedging.py
index 498297440..c9dd53909 100644
--- a/toqito/nonlocal_games/quantum_hedging.py
+++ b/toqito/nonlocal_games/quantum_hedging.py
@@ -48,8 +48,8 @@ class QuantumHedging:
     >>> l_2 = alpha * sin(theta) * e_10
     >>> l_3 = sqrt(1 - alpha ** 2) * cos(theta) * e_01
     >>>
-    >>> q_1 = w_var * w_var.conj().T
-    >>> q_0 = l_1 * l_1.conj().T + l_2 * l_2.conj().T + l_3 * l_3.conj().T
+    >>> q_1 = w_var @ w_var.conj().T
+    >>> q_0 = l_1 @ l_1.conj().T + l_2 @ l_2.conj().T + l_3 @ l_3.conj().T
     >>> molina_watrous = QuantumHedging(q_0, 1)
     >>>
     >>> # cos(pi/8)**2 \approx 0.8536
diff --git a/toqito/perms/permute_systems.py b/toqito/perms/permute_systems.py
index 02ff186df..8cb2c649b 100644
--- a/toqito/perms/permute_systems.py
+++ b/toqito/perms/permute_systems.py
@@ -210,7 +210,7 @@ def permute_systems(
     row_perm = permute_systems(vec_arg, perm, dim[0][:], False, inv_perm)
 
     # This condition is only necessary if the `input_mat` variable is sparse.
-    if isinstance(input_mat, (sparse.csr_matrix, sparse.dia_matrix)):
+    if sparse.issparse(input_mat):
         input_mat = input_mat.toarray()
         permuted_mat = input_mat[row_perm, :]
         permuted_mat = np.array(permuted_mat)
diff --git a/toqito/rand/random_state_vector.py b/toqito/rand/random_state_vector.py
index 69afe8e6f..a7837fb41 100644
--- a/toqito/rand/random_state_vector.py
+++ b/toqito/rand/random_state_vector.py
@@ -26,7 +26,7 @@ def random_state_vector(dim: list[int] | int, is_real: bool = False, k_param: in
     and checking if the density matrix is pure.
 
     >>> from toqito.state_props import is_pure
-    >>> dm = vec.conj().T * vec
+    >>> dm = vec.conj().T @ vec
     >>> is_pure(dm)
     True
 
diff --git a/toqito/state_metrics/helstrom_holevo.py b/toqito/state_metrics/helstrom_holevo.py
index eba976381..868b40aa8 100644
--- a/toqito/state_metrics/helstrom_holevo.py
+++ b/toqito/state_metrics/helstrom_holevo.py
@@ -42,7 +42,7 @@ def helstrom_holevo(rho: np.ndarray, sigma: np.ndarray) -> float:
     >>> e_11 = np.kron(e_1, e_1)
     >>>
     >>> u_vec = 1 / np.sqrt(2) * (e_00 + e_11)
-    >>> rho = u_vec * u_vec.conj().T
+    >>> rho = u_vec @ u_vec.conj().T
     >>> sigma = rho
     >>>
     >>> helstrom_holevo(rho, sigma)
diff --git a/toqito/state_metrics/hilbert_schmidt.py b/toqito/state_metrics/hilbert_schmidt.py
index b119d5bf5..6f7a708b6 100644
--- a/toqito/state_metrics/hilbert_schmidt.py
+++ b/toqito/state_metrics/hilbert_schmidt.py
@@ -23,8 +23,8 @@ def hilbert_schmidt(rho: np.ndarray, sigma: np.ndarray) -> float:
     >>> import numpy as np
     >>> from toqito.states import bell
     >>> from toqito.state_metrics import hilbert_schmidt
-    >>> rho = bell(0) * bell(0).conj().T
-    >>> sigma = bell(3) * bell(3).conj().T
+    >>> rho = bell(0) @ bell(0).conj().T
+    >>> sigma = bell(3) @ bell(3).conj().T
     >>> np.around(hilbert_schmidt(rho, sigma), decimals=2)
     np.float64(1.0)
 
diff --git a/toqito/state_metrics/sub_fidelity.py b/toqito/state_metrics/sub_fidelity.py
index 057b8a26a..b6484d8a1 100644
--- a/toqito/state_metrics/sub_fidelity.py
+++ b/toqito/state_metrics/sub_fidelity.py
@@ -35,8 +35,8 @@ def sub_fidelity(rho: np.ndarray, sigma: np.ndarray) -> float:
     >>> from toqito.states import basis
     >>> from toqito.state_metrics import fidelity
     >>> e_0, e_1 = basis(2, 0), basis(2, 1)
-    >>> rho = 3 / 4 * e_0 * e_0.conj().T + 1 / 4 * e_1 * e_1.conj().T
-    >>> sigma = 1/8 * e_0 * e_0.conj().T + 7/8 * e_1 * e_1.conj().T
+    >>> rho = 3 / 4 * e_0 @ e_0.conj().T + 1 / 4 * e_1 @ e_1.conj().T
+    >>> sigma = 1/8 * e_0 @ e_0.conj().T + 7/8 * e_1 @ e_1.conj().T
     >>> fidelity(rho, sigma)
     np.float64(0.77389339119464)
 
@@ -46,8 +46,8 @@ def sub_fidelity(rho: np.ndarray, sigma: np.ndarray) -> float:
     >>> from toqito.states import basis
     >>> from toqito.state_metrics import sub_fidelity
     >>> e_0, e_1 = basis(2, 0), basis(2, 1)
-    >>> rho = 3 / 4 * e_0 * e_0.conj().T + 1 / 4 * e_1 * e_1.conj().T
-    >>> sigma = 1/8 * e_0 * e_0.conj().T + 7/8 * e_1 * e_1.conj().T
+    >>> rho = 3 / 4 * e_0 @ e_0.conj().T + 1 / 4 * e_1 @ e_1.conj().T
+    >>> sigma = 1/8 * e_0 @ e_0.conj().T + 7/8 * e_1 @ e_1.conj().T
     >>> sub_fidelity(rho, sigma)
     np.float64(0.5989109809347399)
 
@@ -70,5 +70,5 @@ def sub_fidelity(rho: np.ndarray, sigma: np.ndarray) -> float:
         raise ValueError("Sub-fidelity is only defined for density operators.")
 
     return np.real(
-        np.trace(rho * sigma) + np.sqrt(2 * (np.trace(rho * sigma) ** 2 - np.trace(rho * sigma * rho * sigma)))
+        np.trace(rho @ sigma) + np.sqrt(2 * (np.trace(rho @ sigma) ** 2 - np.trace(rho @ sigma @ rho @ sigma)))
     )
diff --git a/toqito/state_metrics/trace_distance.py b/toqito/state_metrics/trace_distance.py
index d1d293a3a..1649011c4 100644
--- a/toqito/state_metrics/trace_distance.py
+++ b/toqito/state_metrics/trace_distance.py
@@ -39,7 +39,7 @@ def trace_distance(rho: np.ndarray, sigma: np.ndarray) -> float:
 
     >>> from toqito.states import bell
     >>> from toqito.state_metrics import trace_distance
-    >>> rho = bell(0) * bell(0).conj().T
+    >>> rho = bell(0) @ bell(0).conj().T
     >>> sigma = rho
     >>> trace_distance(rho, sigma)
     np.float64(0.0)
diff --git a/toqito/state_ops/pure_to_mixed.py b/toqito/state_ops/pure_to_mixed.py
index 4bfe7c3f7..59d6f7e45 100644
--- a/toqito/state_ops/pure_to_mixed.py
+++ b/toqito/state_ops/pure_to_mixed.py
@@ -40,7 +40,7 @@ def pure_to_mixed(phi: np.ndarray) -> np.ndarray:
 
     >>> from toqito.states import bell
     >>> from toqito.state_ops import pure_to_mixed
-    >>> phi = bell(0) * bell(0).conj().T
+    >>> phi = bell(0) @ bell(0).conj().T
     >>> pure_to_mixed(phi)
     array([[0.5, 0. , 0. , 0.5],
            [0. , 0. , 0. , 0. ],
@@ -61,7 +61,7 @@ def pure_to_mixed(phi: np.ndarray) -> np.ndarray:
 
     # It's a pure state vector.
     if min(row_dim, col_dim) == 1:
-        return phi * phi.conj().T
+        return phi @ phi.conj().T
     # It's a density matrix.
     if row_dim == col_dim:
         return phi
diff --git a/toqito/state_opt/optimal_clone.py b/toqito/state_opt/optimal_clone.py
index aac307dfb..26d0346b6 100644
--- a/toqito/state_opt/optimal_clone.py
+++ b/toqito/state_opt/optimal_clone.py
@@ -106,7 +106,7 @@ def optimal_clone(
     # Q = ∑_{k=1}^N p_k |ψ_k ⊗ ψ_k ⊗ ψ_k> <ψ_k ⊗ ψ_k ⊗ ψ_k|
     q_a = np.zeros((dim, dim))
     for k, state in enumerate(states):
-        q_a += probs[k] * tensor(state, state, state.conj()) * tensor(state, state, state.conj()).conj().T
+        q_a += probs[k] * tensor(state, state, state.conj()) @ tensor(state, state, state.conj()).conj().T
 
     # The system is over:
     # Y_1 ⊗ Z_1 ⊗ X_1, ... , Y_n ⊗ Z_n ⊗ X_n.
diff --git a/toqito/state_opt/ppt_distinguishability.py b/toqito/state_opt/ppt_distinguishability.py
index 35f514447..8d2572df3 100644
--- a/toqito/state_opt/ppt_distinguishability.py
+++ b/toqito/state_opt/ppt_distinguishability.py
@@ -79,10 +79,10 @@ def ppt_distinguishability(
     >>> x_4 = np.kron(psi_3, psi_3)
     >>>
     >>> # YDY density matrices.
-    >>> rho_1 = x_1 * x_1.conj().T
-    >>> rho_2 = x_2 * x_2.conj().T
-    >>> rho_3 = x_3 * x_3.conj().T
-    >>> rho_4 = x_4 * x_4.conj().T
+    >>> rho_1 = x_1 @ x_1.conj().T
+    >>> rho_2 = x_2 @ x_2.conj().T
+    >>> rho_3 = x_3 @ x_3.conj().T
+    >>> rho_4 = x_4 @ x_4.conj().T
     >>>
     >>> states = [rho_1, rho_2, rho_3, rho_4]
     >>> probs = [1 / 4, 1 / 4, 1 / 4, 1 / 4]
diff --git a/toqito/state_opt/symmetric_extension_hierarchy.py b/toqito/state_opt/symmetric_extension_hierarchy.py
index 12c3c02c2..8d537083f 100644
--- a/toqito/state_opt/symmetric_extension_hierarchy.py
+++ b/toqito/state_opt/symmetric_extension_hierarchy.py
@@ -89,14 +89,14 @@ def symmetric_extension_hierarchy(
     >>> # Define the resource state.
     >>> eps = 0.5
     >>> eps_state = np.sqrt((1+eps)/2) * e_00 + np.sqrt((1-eps)/2) * e_11
-    >>> eps_dm = eps_state * eps_state.conj().T
+    >>> eps_dm = eps_state @ eps_state.conj().T
     >>>
     >>> # Define the ensemble of states to be distinguished.
     >>> states = [
-    ...     np.kron(bell(0) * bell(0).conj().T, eps_dm),
-    ...     np.kron(bell(1) * bell(1).conj().T, eps_dm),
-    ...     np.kron(bell(2) * bell(2).conj().T, eps_dm),
-    ...     np.kron(bell(3) * bell(3).conj().T, eps_dm),
+    ...     np.kron(bell(0) @ bell(0).conj().T, eps_dm),
+    ...     np.kron(bell(1) @ bell(1).conj().T, eps_dm),
+    ...     np.kron(bell(2) @ bell(2).conj().T, eps_dm),
+    ...     np.kron(bell(3) @ bell(3).conj().T, eps_dm),
     ... ]
     >>>
     >>> # Ensure the distinguishability is conducted on the proper spaces.
@@ -154,7 +154,7 @@ def symmetric_extension_hierarchy(
     # into density matrices.
     if n_cols == 1:
         for i, state_ket in enumerate(states):
-            states[i] = state_ket * state_ket.conj().T
+            states[i] = state_ket @ state_ket.conj().T
 
     # Set default dimension if none was provided.
     if dim is None:
diff --git a/toqito/state_props/entanglement_of_formation.py b/toqito/state_props/entanglement_of_formation.py
index 91e8d0f4c..a6951db12 100644
--- a/toqito/state_props/entanglement_of_formation.py
+++ b/toqito/state_props/entanglement_of_formation.py
@@ -81,7 +81,7 @@ def entanglement_of_formation(rho: np.ndarray, dim: list[int] | int = None) -> f
     if min(dim_x, dim_y) == 1:
         rho = rho[:]
         dim = [int(x) for x in dim]
-        return von_neumann_entropy(partial_trace(rho * rho.conj().T, [1], dim))
+        return von_neumann_entropy(partial_trace(rho @ rho.conj().T, [1], dim))
 
     # Case: :code:`rho` is a density matrix.
     if dim_x == dim_y:
diff --git a/toqito/state_props/has_symmetric_extension.py b/toqito/state_props/has_symmetric_extension.py
index 417430537..da2f2fe6c 100644
--- a/toqito/state_props/has_symmetric_extension.py
+++ b/toqito/state_props/has_symmetric_extension.py
@@ -75,7 +75,7 @@ def has_symmetric_extension(
     >>> from toqito.states import bell
     >>> from toqito.state_props import has_symmetric_extension
     >>>
-    >>> rho = bell(0) * bell(0).conj().T
+    >>> rho = bell(0) @ bell(0).conj().T
     >>> sigma = np.kron(rho, rho)
     >>> has_symmetric_extension(sigma)
     False
diff --git a/toqito/state_props/is_mixed.py b/toqito/state_props/is_mixed.py
index 74ae3d07b..22399653c 100644
--- a/toqito/state_props/is_mixed.py
+++ b/toqito/state_props/is_mixed.py
@@ -27,7 +27,7 @@ def is_mixed(state: np.ndarray) -> bool:
     >>> from toqito.states import basis
     >>> from toqito.state_props import is_mixed
     >>> e_0, e_1 = basis(2, 0), basis(2, 1)
-    >>> rho = 3 / 4 * e_0 * e_0.conj().T + 1 / 4 * e_1 * e_1.conj().T
+    >>> rho = 3 / 4 * e_0 @ e_0.conj().T + 1 / 4 * e_1 @ e_1.conj().T
     >>> is_mixed(rho)
     True
 
diff --git a/toqito/state_props/is_ppt.py b/toqito/state_props/is_ppt.py
index 04ee1b9b1..559b4a101 100644
--- a/toqito/state_props/is_ppt.py
+++ b/toqito/state_props/is_ppt.py
@@ -57,7 +57,7 @@ def is_ppt(mat: np.ndarray, sys: int = 2, dim: int | list[int] = None, tol: floa
 
     >>> from toqito.states import bell
     >>> from toqito.state_props import is_ppt
-    >>> rho = bell(2) * bell(2).conj().T
+    >>> rho = bell(2) @ bell(2).conj().T
     >>> is_ppt(rho)
     False
 
diff --git a/toqito/state_props/is_product.py b/toqito/state_props/is_product.py
index 9ba1775ec..fa12602c8 100644
--- a/toqito/state_props/is_product.py
+++ b/toqito/state_props/is_product.py
@@ -34,7 +34,7 @@ def is_product(rho: np.ndarray, dim: int | list[int] | np.ndarray = None) -> boo
 
     >>> from toqito.state_props import is_product
     >>> from toqito.states import bell
-    >>> rho = bell(0) * bell(0).conj().T
+    >>> rho = bell(0) @ bell(0).conj().T
     >>> u_vec = bell(0)
     >>> is_product(rho)
     (array([False]), None)
diff --git a/toqito/state_props/is_pure.py b/toqito/state_props/is_pure.py
index 36c8eab0f..a6e6d83b8 100644
--- a/toqito/state_props/is_pure.py
+++ b/toqito/state_props/is_pure.py
@@ -36,7 +36,7 @@ def is_pure(state: list[np.ndarray] | np.ndarray) -> bool:
     >>> from toqito.states import bell
     >>> from toqito.state_props import is_pure
     >>> u = bell(0)
-    >>> rho = u * u.conj().T
+    >>> rho = u @ u.conj().T
     >>> is_pure(rho)
     True
 
@@ -47,10 +47,10 @@ def is_pure(state: list[np.ndarray] | np.ndarray) -> bool:
     >>> from toqito.states import bell
     >>> from toqito.state_props import is_pure
     >>> u0, u1, u2, u3 = bell(0), bell(1), bell(2), bell(3)
-    >>> rho0 = u0 * u0.conj().T
-    >>> rho1 = u1 * u1.conj().T
-    >>> rho2 = u2 * u2.conj().T
-    >>> rho3 = u3 * u3.conj().T
+    >>> rho0 = u0 @ u0.conj().T
+    >>> rho1 = u1 @ u1.conj().T
+    >>> rho2 = u2 @ u2.conj().T
+    >>> rho3 = u3 @ u3.conj().T
     >>>
     >>> is_pure([rho0, rho1, rho2, rho3])
     True
diff --git a/toqito/state_props/log_negativity.py b/toqito/state_props/log_negativity.py
index 869bb6a58..73aa244e7 100644
--- a/toqito/state_props/log_negativity.py
+++ b/toqito/state_props/log_negativity.py
@@ -25,7 +25,7 @@ def log_negativity(rho: np.ndarray, dim: list[int] | int = None) -> float:
 
     >>> from toqito.states import bell
     >>> from toqito.state_props import log_negativity
-    >>> rho = bell(0) * bell(0).conj().T
+    >>> rho = bell(0) @ bell(0).conj().T
     >>> log_negativity(rho)
     np.float64(0.9999999999999997)
 
diff --git a/toqito/state_props/negativity.py b/toqito/state_props/negativity.py
index 5bcaff3ec..b46d1bb0b 100644
--- a/toqito/state_props/negativity.py
+++ b/toqito/state_props/negativity.py
@@ -26,7 +26,7 @@ def negativity(rho: np.ndarray, dim: list[int] | int = None) -> float:
 
     >>> from toqito.states import bell
     >>> from toqito.state_props import negativity
-    >>> rho = bell(0) * bell(0).conj().T
+    >>> rho = bell(0) @ bell(0).conj().T
     >>> negativity(rho)
     np.float64(0.4999999999999998)
 
diff --git a/toqito/states/chessboard.py b/toqito/states/chessboard.py
index 67ccb559f..f6bed323e 100644
--- a/toqito/states/chessboard.py
+++ b/toqito/states/chessboard.py
@@ -57,5 +57,5 @@ def chessboard(mat_params: list[float], s_param: float = None, t_param: float =
     v_2 = np.array([[0, mat_params[0], 0, mat_params[1], 0, mat_params[2], 0, 0, 0]])
     v_3 = np.array([[np.conj(mat_params[5]), 0, 0, 0, -np.conj(mat_params[4]), 0, t_param, 0, 0]])
     v_4 = np.array([[0, np.conj(mat_params[1]), 0, -np.conj(mat_params[0]), 0, 0, 0, mat_params[3], 0]])
-    rho = v_1.conj().T * v_1 + v_2.conj().T * v_2 + v_3.conj().T * v_3 + v_4.conj().T * v_4
+    rho = v_1.conj().T @ v_1 + v_2.conj().T @ v_2 + v_3.conj().T @ v_3 + v_4.conj().T @ v_4
     return rho / np.trace(rho)


Test patch:
diff --git a/toqito/channels/tests/test_partial_transpose.py b/toqito/channels/tests/test_partial_transpose.py
index 5cc83ba50..41cbb4b29 100644
--- a/toqito/channels/tests/test_partial_transpose.py
+++ b/toqito/channels/tests/test_partial_transpose.py
@@ -688,7 +688,7 @@ def test_partial_transpose_16_by_16_subsystems_2_2_2_2():
 
 def test_partial_transpose_bell_state():
     """Test partial transpose on a Bell state."""
-    rho = bell(2) * bell(2).conj().T
+    rho = bell(2) @ bell(2).conj().T
     expected_res = np.array([[0, 0, 0, 1 / 2], [0, 1 / 2, 0, 0], [0, 0, 1 / 2, 0], [1 / 2, 0, 0, 0]])
     res = partial_transpose(rho)
     np.testing.assert_equal(np.allclose(res, expected_res), True)
diff --git a/toqito/matrix_props/tests/test_majorizes.py b/toqito/matrix_props/tests/test_majorizes.py
index 814416b32..7ad21cdde 100644
--- a/toqito/matrix_props/tests/test_majorizes.py
+++ b/toqito/matrix_props/tests/test_majorizes.py
@@ -15,12 +15,12 @@ def test_majorizes_simple_example():
 def test_majorizes_max_entangled():
     """Test that max entangled partial trace returns False."""
     v_vec = max_entangled(3)
-    rho = v_vec * v_vec.conj().T
+    rho = v_vec @ v_vec.conj().T
     np.testing.assert_equal(majorizes(partial_trace(rho, [1], [3, 3]), rho), False)
 
 
 def test_majorizes_max_entangled_flip():
     """Test that max entangled partial trace returns True (flipped args)."""
     v_vec = max_entangled(3)
-    rho = v_vec * v_vec.conj().T
+    rho = v_vec @ v_vec.conj().T
     np.testing.assert_equal(majorizes(rho, partial_trace(rho, [1], [3, 3])), True)
diff --git a/toqito/matrix_props/tests/test_spark.py b/toqito/matrix_props/tests/test_spark.py
new file mode 100644
index 000000000..bf0d7b09f
--- /dev/null
+++ b/toqito/matrix_props/tests/test_spark.py
@@ -0,0 +1,38 @@
+"""Tests for spark function."""
+
+import numpy as np
+import pytest
+
+from toqito.matrix_props import spark
+
+
+@pytest.mark.parametrize("matrix, expected_spark, description", [
+    (np.array([[1, 0, 1], [0, 1, 1], [1, 1, 0]]), 4, "Square matrix with all columns linearly independent"),
+    (np.array([[1, 0, 1, 2], [0, 1, 1, 3]]), 3, "Non-square matrix with first three columns linearly dependent"),
+    (np.array([[1, 0, 0], [1, 1, 0], [1, 1, 0]]), 1, "Matrix with a zero column"),
+    (np.array([[1, 0, 0], [0, 1, 0], [0, 0, 1]]), 4, "Full rank matrix"),
+    (np.array([[1, 2, 3, 4], [0, 1, 1, 2], [1, 0, 1, 1]]), 3, "Matrix with linearly dependent subset of columns"),
+    (np.array([[1, 0, 0, 1, 2],
+               [0, 1, 0, 1, 1],
+               [0, 0, 1, 1, 0],
+               [1, 1, 1, 0, 1]]), 5, "Larger matrix with all columns needed for linear dependence"),
+])
+def test_spark(matrix, expected_spark, description):
+    """Test spark function with various input matrices."""
+    assert spark(matrix) == expected_spark, f"Failed for case: {description}"
+
+def test_spark_property_rank():
+    """Test spark function property: spark(A) <= rank(A) + 1."""
+    A = np.random.rand(3, 5)
+    s = spark(A)
+    r = np.linalg.matrix_rank(A)
+    assert s <= r + 1, "Spark should be <= rank + 1"
+
+@pytest.mark.parametrize("invalid_input", [
+    np.array([1, 2, 3]),  # 1D array
+    np.random.rand(2, 2, 2),  # 3D array
+])
+def test_spark_invalid_input(invalid_input):
+    """Test spark function with invalid inputs."""
+    with pytest.raises(ValueError):
+        spark(invalid_input)
diff --git a/toqito/matrix_props/tests/test_trace_norm.py b/toqito/matrix_props/tests/test_trace_norm.py
index f7e7f7b85..bccd0c9cc 100644
--- a/toqito/matrix_props/tests/test_trace_norm.py
+++ b/toqito/matrix_props/tests/test_trace_norm.py
@@ -13,7 +13,7 @@ def test_trace_norm():
     e_11 = np.kron(e_1, e_1)
 
     u_vec = 1 / np.sqrt(2) * (e_00 + e_11)
-    rho = u_vec * u_vec.conj().T
+    rho = u_vec @ u_vec.conj().T
 
     res = trace_norm(rho)
     _, singular_vals, _ = np.linalg.svd(rho)
diff --git a/toqito/measurement_ops/tests/test_measure.py b/toqito/measurement_ops/tests/test_measure.py
index 457492354..31c59b62f 100644
--- a/toqito/measurement_ops/tests/test_measure.py
+++ b/toqito/measurement_ops/tests/test_measure.py
@@ -10,9 +10,9 @@ def test_measure_state():
     """Test measure on quantum state."""
     e_0, e_1 = basis(2, 0), basis(2, 1)
     psi = 1 / np.sqrt(3) * e_0 + np.sqrt(2 / 3) * e_1
-    rho = psi * psi.conj().T
+    rho = psi @ psi.conj().T
 
-    proj_0 = e_0 * e_0.conj().T
-    proj_1 = e_1 * e_1.conj().T
+    proj_0 = e_0 @ e_0.conj().T
+    proj_1 = e_1 @ e_1.conj().T
     np.testing.assert_equal(np.isclose(measure(proj_0, rho), 1 / 3), True)
     np.testing.assert_equal(np.isclose(measure(proj_1, rho), 2 / 3), True)
diff --git a/toqito/nonlocal_games/tests/test_extended_nonlocal_game.py b/toqito/nonlocal_games/tests/test_extended_nonlocal_game.py
index 768a3e1b2..b0370c505 100644
--- a/toqito/nonlocal_games/tests/test_extended_nonlocal_game.py
+++ b/toqito/nonlocal_games/tests/test_extended_nonlocal_game.py
@@ -23,10 +23,10 @@ def bb84_extended_nonlocal_game():
         num_alice_in, num_bob_in = 2, 2
 
         pred_mat = np.zeros([dim, dim, num_alice_out, num_bob_out, num_alice_in, num_bob_in])
-        pred_mat[:, :, 0, 0, 0, 0] = e_0 * e_0.conj().T
-        pred_mat[:, :, 0, 0, 1, 1] = e_p * e_p.conj().T
-        pred_mat[:, :, 1, 1, 0, 0] = e_1 * e_1.conj().T
-        pred_mat[:, :, 1, 1, 1, 1] = e_m * e_m.conj().T
+        pred_mat[:, :, 0, 0, 0, 0] = e_0 @ e_0.conj().T
+        pred_mat[:, :, 0, 0, 1, 1] = e_p @ e_p.conj().T
+        pred_mat[:, :, 1, 1, 0, 0] = e_1 @ e_1.conj().T
+        pred_mat[:, :, 1, 1, 1, 1] = e_m @ e_m.conj().T
 
         prob_mat = 1 / 2 * np.identity(2)
 
@@ -88,21 +88,21 @@ def moe_mub_4_in_3_out_game():
         num_out = 3
         pred_mat = np.zeros([dim, dim, num_out, num_out, num_in, num_in], dtype=complex)
 
-        pred_mat[:, :, 0, 0, 0, 0] = mubs[0][0] * mubs[0][0].conj().T
-        pred_mat[:, :, 1, 1, 0, 0] = mubs[0][1] * mubs[0][1].conj().T
-        pred_mat[:, :, 2, 2, 0, 0] = mubs[0][2] * mubs[0][2].conj().T
+        pred_mat[:, :, 0, 0, 0, 0] = mubs[0][0] @ mubs[0][0].conj().T
+        pred_mat[:, :, 1, 1, 0, 0] = mubs[0][1] @ mubs[0][1].conj().T
+        pred_mat[:, :, 2, 2, 0, 0] = mubs[0][2] @ mubs[0][2].conj().T
 
-        pred_mat[:, :, 0, 0, 1, 1] = mubs[1][0] * mubs[1][0].conj().T
-        pred_mat[:, :, 1, 1, 1, 1] = mubs[1][1] * mubs[1][1].conj().T
-        pred_mat[:, :, 2, 2, 1, 1] = mubs[1][2] * mubs[1][2].conj().T
+        pred_mat[:, :, 0, 0, 1, 1] = mubs[1][0] @ mubs[1][0].conj().T
+        pred_mat[:, :, 1, 1, 1, 1] = mubs[1][1] @ mubs[1][1].conj().T
+        pred_mat[:, :, 2, 2, 1, 1] = mubs[1][2] @ mubs[1][2].conj().T
 
-        pred_mat[:, :, 0, 0, 2, 2] = mubs[2][0] * mubs[2][0].conj().T
-        pred_mat[:, :, 1, 1, 2, 2] = mubs[2][1] * mubs[2][1].conj().T
-        pred_mat[:, :, 2, 2, 2, 2] = mubs[2][2] * mubs[2][2].conj().T
+        pred_mat[:, :, 0, 0, 2, 2] = mubs[2][0] @ mubs[2][0].conj().T
+        pred_mat[:, :, 1, 1, 2, 2] = mubs[2][1] @ mubs[2][1].conj().T
+        pred_mat[:, :, 2, 2, 2, 2] = mubs[2][2] @ mubs[2][2].conj().T
 
-        pred_mat[:, :, 0, 0, 3, 3] = mubs[3][0] * mubs[3][0].conj().T
-        pred_mat[:, :, 1, 1, 3, 3] = mubs[3][1] * mubs[3][1].conj().T
-        pred_mat[:, :, 2, 2, 3, 3] = mubs[3][2] * mubs[3][2].conj().T
+        pred_mat[:, :, 0, 0, 3, 3] = mubs[3][0] @ mubs[3][0].conj().T
+        pred_mat[:, :, 1, 1, 3, 3] = mubs[3][1] @ mubs[3][1].conj().T
+        pred_mat[:, :, 2, 2, 3, 3] = mubs[3][2] @ mubs[3][2].conj().T
 
         return prob_mat, pred_mat
 
diff --git a/toqito/nonlocal_games/tests/test_quantum_hedging.py b/toqito/nonlocal_games/tests/test_quantum_hedging.py
index 4c2d1a117..63f2af7bd 100644
--- a/toqito/nonlocal_games/tests/test_quantum_hedging.py
+++ b/toqito/nonlocal_games/tests/test_quantum_hedging.py
@@ -26,8 +26,8 @@ class TestQuantumHedging(unittest.TestCase):
 
     l_3 = sqrt(1 - alpha**2) * cos(theta) * e_01
 
-    q_1 = w_var * w_var.conj().T
-    q_0 = l_1 * l_1.conj().T + l_2 * l_2.conj().T + l_3 * l_3.conj().T
+    q_1 = w_var @ w_var.conj().T
+    q_0 = l_1 @ l_1.conj().T + l_2 @ l_2.conj().T + l_3 @ l_3.conj().T
 
     def test_max_prob_outcome_a_primal_1_dim(self):
         """Maximal probability of outcome "a" when dim == 1.
diff --git a/toqito/state_metrics/tests/test_bures_angle.py b/toqito/state_metrics/tests/test_bures_angle.py
index ff998dd64..072c5daa4 100644
--- a/toqito/state_metrics/tests/test_bures_angle.py
+++ b/toqito/state_metrics/tests/test_bures_angle.py
@@ -18,8 +18,8 @@ def test_bures_angle_default():
 def test_bures_angle_non_identical_states_1():
     """Test the bures_angle between two non-identical states."""
     e_0, e_1 = basis(2, 0), basis(2, 1)
-    rho = 3 / 4 * e_0 * e_0.conj().T + 1 / 4 * e_1 * e_1.conj().T
-    sigma = 2 / 3 * e_0 * e_0.conj().T + 1 / 3 * e_1 * e_1.conj().T
+    rho = 3 / 4 * e_0 @ e_0.conj().T + 1 / 4 * e_1 @ e_1.conj().T
+    sigma = 2 / 3 * e_0 @ e_0.conj().T + 1 / 3 * e_1 @ e_1.conj().T
 
     ang = bures_angle(rho, sigma)
     np.testing.assert_equal(np.isclose(ang, 0.06499, rtol=1e-03), True)
@@ -28,8 +28,8 @@ def test_bures_angle_non_identical_states_1():
 def test_bures_angle_non_identical_states_2():
     """Test the bures_angle between two non-identical states."""
     e_0, e_1 = basis(2, 0), basis(2, 1)
-    rho = 3 / 4 * e_0 * e_0.conj().T + 1 / 4 * e_1 * e_1.conj().T
-    sigma = 1 / 8 * e_0 * e_0.conj().T + 7 / 8 * e_1 * e_1.conj().T
+    rho = 3 / 4 * e_0 @ e_0.conj().T + 1 / 4 * e_1 @ e_1.conj().T
+    sigma = 1 / 8 * e_0 @ e_0.conj().T + 7 / 8 * e_1 @ e_1.conj().T
 
     ang = bures_angle(rho, sigma)
     np.testing.assert_equal(np.isclose(ang, 0.4955, rtol=1e-03), True)
@@ -39,8 +39,8 @@ def test_bures_angle_pure_states():
     """Test the bures_angle between two pure states."""
     e_0, e_1 = basis(2, 0), basis(2, 1)
     e_plus = (e_0 + e_1) / np.sqrt(2)
-    rho = e_plus * e_plus.conj().T
-    sigma = e_0 * e_0.conj().T
+    rho = e_plus @ e_plus.conj().T
+    sigma = e_0 @ e_0.conj().T
 
     ang = bures_angle(rho, sigma)
     np.testing.assert_equal(np.isclose(ang, 0.5718, rtol=1e-03), True)
diff --git a/toqito/state_metrics/tests/test_bures_distance.py b/toqito/state_metrics/tests/test_bures_distance.py
index 3808e0a9c..c27036bee 100644
--- a/toqito/state_metrics/tests/test_bures_distance.py
+++ b/toqito/state_metrics/tests/test_bures_distance.py
@@ -18,16 +18,16 @@ def test_bures_distance_default():
 def test_bures_distance_non_identical_states_1():
     """Test the bures_distance between two non-identical states."""
     e_0, e_1 = basis(2, 0), basis(2, 1)
-    rho = 3 / 4 * e_0 * e_0.conj().T + 1 / 4 * e_1 * e_1.conj().T
-    sigma = 2 / 3 * e_0 * e_0.conj().T + 1 / 3 * e_1 * e_1.conj().T
+    rho = 3 / 4 * e_0 @ e_0.conj().T + 1 / 4 * e_1 @ e_1.conj().T
+    sigma = 2 / 3 * e_0 @ e_0.conj().T + 1 / 3 * e_1 @ e_1.conj().T
     np.testing.assert_equal(np.isclose(bures_distance(rho, sigma), 0.0918, rtol=1e-03), True)
 
 
 def test_bures_distance_non_identical_states_2():
     """Test the bures_distance between two non-identical states."""
     e_0, e_1 = basis(2, 0), basis(2, 1)
-    rho = 3 / 4 * e_0 * e_0.conj().T + 1 / 4 * e_1 * e_1.conj().T
-    sigma = 1 / 8 * e_0 * e_0.conj().T + 7 / 8 * e_1 * e_1.conj().T
+    rho = 3 / 4 * e_0 @ e_0.conj().T + 1 / 4 * e_1 @ e_1.conj().T
+    sigma = 1 / 8 * e_0 @ e_0.conj().T + 7 / 8 * e_1 @ e_1.conj().T
     np.testing.assert_equal(np.isclose(bures_distance(rho, sigma), 0.6724, rtol=1e-03), True)
 
 
@@ -35,8 +35,8 @@ def test_bures_distance_pure_states():
     """Test the bures_distance between two pure states."""
     e_0, e_1 = basis(2, 0), basis(2, 1)
     e_plus = (e_0 + e_1) / np.sqrt(2)
-    rho = e_plus * e_plus.conj().T
-    sigma = e_0 * e_0.conj().T
+    rho = e_plus @ e_plus.conj().T
+    sigma = e_0 @ e_0.conj().T
     np.testing.assert_equal(np.isclose(bures_distance(rho, sigma), 0.765, rtol=1e-03), True)
 
 
diff --git a/toqito/state_metrics/tests/test_fidelity.py b/toqito/state_metrics/tests/test_fidelity.py
index 2fd4aa82b..1b6e4394a 100644
--- a/toqito/state_metrics/tests/test_fidelity.py
+++ b/toqito/state_metrics/tests/test_fidelity.py
@@ -28,16 +28,16 @@ def test_fidelity_cvx():
 def test_fidelity_non_identical_states_1():
     """Test the fidelity between two non-identical states."""
     e_0, e_1 = basis(2, 0), basis(2, 1)
-    rho = 3 / 4 * e_0 * e_0.conj().T + 1 / 4 * e_1 * e_1.conj().T
-    sigma = 2 / 3 * e_0 * e_0.conj().T + 1 / 3 * e_1 * e_1.conj().T
+    rho = 3 / 4 * e_0 @ e_0.conj().T + 1 / 4 * e_1 @ e_1.conj().T
+    sigma = 2 / 3 * e_0 @ e_0.conj().T + 1 / 3 * e_1 @ e_1.conj().T
     np.testing.assert_equal(np.isclose(fidelity(rho, sigma), 0.996, rtol=1e-03), True)
 
 
 def test_fidelity_non_identical_states_2():
     """Test the fidelity between two non-identical states."""
     e_0, e_1 = basis(2, 0), basis(2, 1)
-    rho = 3 / 4 * e_0 * e_0.conj().T + 1 / 4 * e_1 * e_1.conj().T
-    sigma = 1 / 8 * e_0 * e_0.conj().T + 7 / 8 * e_1 * e_1.conj().T
+    rho = 3 / 4 * e_0 @ e_0.conj().T + 1 / 4 * e_1 @ e_1.conj().T
+    sigma = 1 / 8 * e_0 @ e_0.conj().T + 7 / 8 * e_1 @ e_1.conj().T
     np.testing.assert_equal(np.isclose(fidelity(rho, sigma), 0.774, rtol=1e-03), True)
 
 
diff --git a/toqito/state_metrics/tests/test_helstrom_holevo.py b/toqito/state_metrics/tests/test_helstrom_holevo.py
index d6054eacf..a464262f2 100644
--- a/toqito/state_metrics/tests/test_helstrom_holevo.py
+++ b/toqito/state_metrics/tests/test_helstrom_holevo.py
@@ -13,7 +13,7 @@ def test_helstrom_holevo_same_state():
     e_11 = np.kron(e_1, e_1)
 
     u_vec = 1 / np.sqrt(2) * (e_00 + e_11)
-    sigma_0 = u_vec * u_vec.conj().T
+    sigma_0 = u_vec @ u_vec.conj().T
     sigma_1 = sigma_0
 
     res = helstrom_holevo(sigma_0, sigma_1)
diff --git a/toqito/state_metrics/tests/test_hilbert_schmidt.py b/toqito/state_metrics/tests/test_hilbert_schmidt.py
index 8aca89593..6b5aa8448 100644
--- a/toqito/state_metrics/tests/test_hilbert_schmidt.py
+++ b/toqito/state_metrics/tests/test_hilbert_schmidt.py
@@ -8,8 +8,8 @@
 
 def test_hilbert_schmidt_bell():
     r"""Test Hilbert-Schmidt distance on two Bell states."""
-    rho = bell(0) * bell(0).conj().T
-    sigma = bell(3) * bell(3).conj().T
+    rho = bell(0) @ bell(0).conj().T
+    sigma = bell(3) @ bell(3).conj().T
 
     res = hilbert_schmidt(rho, sigma)
 
diff --git a/toqito/state_metrics/tests/test_matsumoto_fidelity.py b/toqito/state_metrics/tests/test_matsumoto_fidelity.py
index f89a4b686..0f9367986 100644
--- a/toqito/state_metrics/tests/test_matsumoto_fidelity.py
+++ b/toqito/state_metrics/tests/test_matsumoto_fidelity.py
@@ -10,11 +10,11 @@
 rho = np.array([[1 / 2, 0, 0, 1 / 2], [0, 0, 0, 0], [0, 0, 0, 0], [1 / 2, 0, 0, 1 / 2]])
 
 e_0, e_1 = basis(2, 0), basis(2, 1)
-rho1 = 3 / 4 * e_0 * e_0.conj().T + 1 / 4 * e_1 * e_1.conj().T
-sigma1 = 2 / 3 * e_0 * e_0.conj().T + 1 / 3 * e_1 * e_1.conj().T
+rho1 = 3 / 4 * e_0 @ e_0.conj().T + 1 / 4 * e_1 @ e_1.conj().T
+sigma1 = 2 / 3 * e_0 @ e_0.conj().T + 1 / 3 * e_1 @ e_1.conj().T
 
-rho2 = 3 / 4 * e_0 * e_0.conj().T + 1 / 4 * e_1 * e_1.conj().T
-sigma2 = 1 / 8 * e_0 * e_0.conj().T + 7 / 8 * e_1 * e_1.conj().T
+rho2 = 3 / 4 * e_0 @ e_0.conj().T + 1 / 4 * e_1 @ e_1.conj().T
+sigma2 = 1 / 8 * e_0 @ e_0.conj().T + 7 / 8 * e_1 @ e_1.conj().T
 
 rho5 = cvxpy.bmat([[1 / 2, 0, 0, 1 / 2], [0, 0, 0, 0], [0, 0, 0, 0], [1 / 2, 0, 0, 1 / 2]])
 
diff --git a/toqito/state_metrics/tests/test_sub_fidelity.py b/toqito/state_metrics/tests/test_sub_fidelity.py
index 35dc6943b..d14be1f6a 100644
--- a/toqito/state_metrics/tests/test_sub_fidelity.py
+++ b/toqito/state_metrics/tests/test_sub_fidelity.py
@@ -8,7 +8,7 @@
 
 def test_sub_fidelity_default():
     """Test sub_fidelity default arguments."""
-    rho = bell(0) * bell(0).conj().T
+    rho = bell(0) @ bell(0).conj().T
     sigma = rho
 
     res = sub_fidelity(rho, sigma)
@@ -18,8 +18,8 @@ def test_sub_fidelity_default():
 def test_sub_fidelity_lower_bound_1():
     """Test sub_fidelity is lower bound on fidelity for rho and sigma."""
     e_0, e_1 = basis(2, 0), basis(2, 1)
-    rho = 3 / 4 * e_0 * e_0.conj().T + 1 / 4 * e_1 * e_1.conj().T
-    sigma = 2 / 3 * e_0 * e_0.conj().T + 1 / 3 * e_1 * e_1.conj().T
+    rho = 3 / 4 * e_0 @ e_0.conj().T + 1 / 4 * e_1 @ e_1.conj().T
+    sigma = 2 / 3 * e_0 @ e_0.conj().T + 1 / 3 * e_1 @ e_1.conj().T
 
     res = sub_fidelity(rho, sigma)
     np.testing.assert_array_less(res, fidelity(rho, sigma))
@@ -28,8 +28,8 @@ def test_sub_fidelity_lower_bound_1():
 def test_sub_fidelity_lower_bound_2():
     """Test sub_fidelity is lower bound on fidelity for rho and pi."""
     e_0, e_1 = basis(2, 0), basis(2, 1)
-    rho = 3 / 4 * e_0 * e_0.conj().T + 1 / 4 * e_1 * e_1.conj().T
-    sigma = 1 / 8 * e_0 * e_0.conj().T + 7 / 8 * e_1 * e_1.conj().T
+    rho = 3 / 4 * e_0 @ e_0.conj().T + 1 / 4 * e_1 @ e_1.conj().T
+    sigma = 1 / 8 * e_0 @ e_0.conj().T + 7 / 8 * e_1 @ e_1.conj().T
 
     res = sub_fidelity(rho, sigma)
     np.testing.assert_array_less(res, fidelity(rho, sigma))
diff --git a/toqito/state_metrics/tests/test_trace_distance.py b/toqito/state_metrics/tests/test_trace_distance.py
index d2ddf1169..ec20c458d 100644
--- a/toqito/state_metrics/tests/test_trace_distance.py
+++ b/toqito/state_metrics/tests/test_trace_distance.py
@@ -13,7 +13,7 @@ def test_trace_distance_same_state():
     e_11 = np.kron(e_1, e_1)
 
     u_vec = 1 / np.sqrt(2) * (e_00 + e_11)
-    rho = u_vec * u_vec.conj().T
+    rho = u_vec @ u_vec.conj().T
     sigma = rho
 
     res = trace_distance(rho, sigma)
diff --git a/toqito/state_ops/tests/test_pure_to_mixed.py b/toqito/state_ops/tests/test_pure_to_mixed.py
index 187b7df27..63784a596 100644
--- a/toqito/state_ops/tests/test_pure_to_mixed.py
+++ b/toqito/state_ops/tests/test_pure_to_mixed.py
@@ -21,7 +21,7 @@ def test_pure_to_mixed_density_matrix():
     """Convert pure state to mixed state density matrix."""
     expected_res = np.array([[1 / 2, 0, 0, 1 / 2], [0, 0, 0, 0], [0, 0, 0, 0], [1 / 2, 0, 0, 1 / 2]])
 
-    phi = bell(0) * bell(0).conj().T
+    phi = bell(0) @ bell(0).conj().T
     res = pure_to_mixed(phi)
 
     bool_mat = np.isclose(res, expected_res)
diff --git a/toqito/state_opt/tests/test_ppt_distinguishability.py b/toqito/state_opt/tests/test_ppt_distinguishability.py
index 39f8fff87..76e34759c 100644
--- a/toqito/state_opt/tests/test_ppt_distinguishability.py
+++ b/toqito/state_opt/tests/test_ppt_distinguishability.py
@@ -22,10 +22,10 @@ def test_ppt_distinguishability_yyd_density_matrices():
     x_3 = np.kron(psi_2, psi_3)
     x_4 = np.kron(psi_3, psi_3)
 
-    rho_1 = x_1 * x_1.conj().T
-    rho_2 = x_2 * x_2.conj().T
-    rho_3 = x_3 * x_3.conj().T
-    rho_4 = x_4 * x_4.conj().T
+    rho_1 = x_1 @ x_1.conj().T
+    rho_2 = x_2 @ x_2.conj().T
+    rho_3 = x_3 @ x_3.conj().T
+    rho_4 = x_4 @ x_4.conj().T
 
     states = [rho_1, rho_2, rho_3, rho_4]
     probs = [1 / 4, 1 / 4, 1 / 4, 1 / 4]
@@ -128,10 +128,10 @@ def test_ppt_distinguishability_yyd_states_no_probs():
     x_3 = np.kron(psi_2, psi_3)
     x_4 = np.kron(psi_3, psi_3)
 
-    rho_1 = x_1 * x_1.conj().T
-    rho_2 = x_2 * x_2.conj().T
-    rho_3 = x_3 * x_3.conj().T
-    rho_4 = x_4 * x_4.conj().T
+    rho_1 = x_1 @ x_1.conj().T
+    rho_2 = x_2 @ x_2.conj().T
+    rho_3 = x_3 @ x_3.conj().T
+    rho_4 = x_4 @ x_4.conj().T
 
     states = [rho_1, rho_2, rho_3, rho_4]
 
diff --git a/toqito/state_opt/tests/test_symmetric_extension_hierarchy.py b/toqito/state_opt/tests/test_symmetric_extension_hierarchy.py
index 62807ded6..2ad3465d8 100644
--- a/toqito/state_opt/tests/test_symmetric_extension_hierarchy.py
+++ b/toqito/state_opt/tests/test_symmetric_extension_hierarchy.py
@@ -11,10 +11,10 @@
 def test_symmetric_extension_hierarchy_four_bell_density_matrices():
     """Symmetric extension hierarchy for four Bell density matrices."""
     states = [
-        bell(0) * bell(0).conj().T,
-        bell(1) * bell(1).conj().T,
-        bell(2) * bell(2).conj().T,
-        bell(3) * bell(3).conj().T,
+        bell(0) @ bell(0).conj().T,
+        bell(1) @ bell(1).conj().T,
+        bell(2) @ bell(2).conj().T,
+        bell(3) @ bell(3).conj().T,
     ]
     res = symmetric_extension_hierarchy(states=states, probs=None, level=2)
     np.testing.assert_equal(np.isclose(res, 1 / 2, atol=1e-5), True)
@@ -37,10 +37,10 @@ def test_symmetric_extension_hierarchy_four_bell_with_resource_state_lvl_1():
     eps_dm = eps_state * eps_state.conj().T
 
     states = [
-        np.kron(bell(0) * bell(0).conj().T, eps_dm),
-        np.kron(bell(1) * bell(1).conj().T, eps_dm),
-        np.kron(bell(2) * bell(2).conj().T, eps_dm),
-        np.kron(bell(3) * bell(3).conj().T, eps_dm),
+        np.kron(bell(0) @ bell(0).conj().T, eps_dm),
+        np.kron(bell(1) @ bell(1).conj().T, eps_dm),
+        np.kron(bell(2) @ bell(2).conj().T, eps_dm),
+        np.kron(bell(3) @ bell(3).conj().T, eps_dm),
     ]
 
     # Ensure we are checking the correct partition of the states.
@@ -70,10 +70,10 @@ def test_symmetric_extension_hierarchy_four_bell_with_resource_state():
     eps_dm = eps_state * eps_state.conj().T
 
     states = [
-        np.kron(bell(0) * bell(0).conj().T, eps_dm),
-        np.kron(bell(1) * bell(1).conj().T, eps_dm),
-        np.kron(bell(2) * bell(2).conj().T, eps_dm),
-        np.kron(bell(3) * bell(3).conj().T, eps_dm),
+        np.kron(bell(0) @ bell(0).conj().T, eps_dm),
+        np.kron(bell(1) @ bell(1).conj().T, eps_dm),
+        np.kron(bell(2) @ bell(2).conj().T, eps_dm),
+        np.kron(bell(3) @ bell(3).conj().T, eps_dm),
     ]
 
     # Ensure we are checking the correct partition of the states.
@@ -98,12 +98,12 @@ def test_symmetric_extension_hierarchy_three_bell_with_resource_state():
 
     eps = 1 / 2
     eps_state = np.sqrt((1 + eps) / 2) * e_00 + np.sqrt((1 - eps) / 2) * e_11
-    eps_dm = eps_state * eps_state.conj().T
+    eps_dm = eps_state @ eps_state.conj().T
 
     states = [
-        np.kron(bell(0) * bell(0).conj().T, eps_dm),
-        np.kron(bell(1) * bell(1).conj().T, eps_dm),
-        np.kron(bell(2) * bell(2).conj().T, eps_dm),
+        np.kron(bell(0) @ bell(0).conj().T, eps_dm),
+        np.kron(bell(1) @ bell(1).conj().T, eps_dm),
+        np.kron(bell(2) @ bell(2).conj().T, eps_dm),
     ]
 
     # Ensure we are checking the correct partition of the states.
@@ -120,8 +120,8 @@ def test_symmetric_extension_hierarchy_three_bell_with_resource_state():
 def test_invalid_symmetric_extension_hierarchy_probs():
     """Invalid probability vector for symmetric extension hierarchy."""
     with np.testing.assert_raises(ValueError):
-        rho1 = bell(0) * bell(0).conj().T
-        rho2 = bell(1) * bell(1).conj().T
+        rho1 = bell(0) @ bell(0).conj().T
+        rho2 = bell(1) @ bell(1).conj().T
         states = [rho1, rho2]
         symmetric_extension_hierarchy(states, [1, 2, 3])
 
diff --git a/toqito/state_props/tests/test_is_separable.py b/toqito/state_props/tests/test_is_separable.py
index 14330ae32..6c76f153f 100644
--- a/toqito/state_props/tests/test_is_separable.py
+++ b/toqito/state_props/tests/test_is_separable.py
@@ -61,7 +61,7 @@ def test_invalid_dim_parameter():
 
 def test_entangled_ppt_criterion():
     """Determined to be entangled via the PPT criterion."""
-    rho = bell(0) * bell(0).conj().T
+    rho = bell(0) @ bell(0).conj().T
     np.testing.assert_equal(is_separable(rho), False)
 
 
@@ -72,7 +72,7 @@ def test_ppt_small_dimensions():
 
     e_0, e_1 = basis(2, 0), basis(2, 1)
     phi = np.kron((1 / np.sqrt(2) * e_0 + 1 / np.sqrt(2) * e_1), psi)
-    sigma = phi * phi.conj().T
+    sigma = phi @ phi.conj().T
     np.testing.assert_equal(is_separable(sigma), True)
 
 
@@ -102,7 +102,7 @@ def test_entangled_realignment_criterion():
     # :math:`\rho = \frac{1}{4} \mathbb{I}_3 \otimes \mathbb{I}_3 - \sum_{i=0}^4 | \psi_i \rangle \langle \psi_i |`
     rho = np.identity(9)
     for i in range(5):
-        rho = rho - tile(i) * tile(i).conj().T
+        rho = rho - tile(i) @ tile(i).conj().T
     rho = rho / 4
     np.testing.assert_equal(is_density(rho), True)
     np.testing.assert_equal(is_separable(rho), False)
diff --git a/toqito/states/tests/test_singlet.py b/toqito/states/tests/test_singlet.py
index dc34bda6a..0bdf0e81b 100644
--- a/toqito/states/tests/test_singlet.py
+++ b/toqito/states/tests/test_singlet.py
@@ -9,7 +9,7 @@ def test_gen_bell_dim_2():
     """Generalized singlet state for dim = 2."""
     dim = 2
 
-    expected_res = bell(3) * bell(3).conj().T
+    expected_res = bell(3) @ bell(3).conj().T
 
     res = singlet(dim)
     np.testing.assert_allclose(res, expected_res)