Repository: vprusso/toqito

Issue:
Feature: Generate random orthonormal basis
This particular feature may be of use if we wish to generate a collection of quantum states (basis of vectors) that have the property of being orthonormal. The columns of a randomly generated unitary matrix can be used to generate an orthonormal basis

Something like this could be used:

```py
def random_orthonormal_basis(dim: int) -> list[np.ndarray]:
    # Generate a random unitary matrix
    U = unitary_group.rvs(dim)

    # The columns of U form an orthonormal basis
    basis = [U[:, i] for i in range(dim)]
    return basis
```

This feature should live in `rand/random_orthonormal_basis.py` and would need to have the requisite coverage and testing in `tests/test_random_orthonormal_basis.py`.

Hints from issue comments:

Solution patch:
diff --git a/docs/refs.bib b/docs/refs.bib
index 5df631a7c..2921eafa9 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,
@@ -806,6 +813,15 @@ @misc{Russo_2017_Extended
 }
 
 #Last name begins with S
+
+@misc{SE_1688950,
+  author  = "Stack Exchange Mathematics",
+  title   = "Why do the columns of a unitary matrix form an orthonormal basis?",
+  howpublished = {https://math.stackexchange.com/q/1688950}
+  
+}
+
+
 @misc{Seshadri_2021_Git,
   author  = "Seshadri, Akshay",
   title   = "Minimax Fidelity Estimation",
diff --git a/poetry.lock b/poetry.lock
index 53ae6a79f..a8ae5bb53 100644
--- a/poetry.lock
+++ b/poetry.lock
@@ -611,13 +611,13 @@ files = [
 
 [[package]]
 name = "ipython"
-version = "8.27.0"
+version = "8.28.0"
 description = "IPython: Productive Interactive Computing"
 optional = false
 python-versions = ">=3.10"
 files = [
-    {file = "ipython-8.27.0-py3-none-any.whl", hash = "sha256:f68b3cb8bde357a5d7adc9598d57e22a45dfbea19eb6b98286fa3b288c9cd55c"},
-    {file = "ipython-8.27.0.tar.gz", hash = "sha256:0b99a2dc9f15fd68692e898e5568725c6d49c527d36a9fb5960ffbdeaa82ff7e"},
+    {file = "ipython-8.28.0-py3-none-any.whl", hash = "sha256:530ef1e7bb693724d3cdc37287c80b07ad9b25986c007a53aa1857272dac3f35"},
+    {file = "ipython-8.28.0.tar.gz", hash = "sha256:0d0d15ca1e01faeb868ef56bc7ee5a0de5bd66885735682e8a322ae289a13d1a"},
 ]
 
 [package.dependencies]
@@ -2034,13 +2034,13 @@ dev = ["hypothesis (>=6.70.0)", "pytest (>=7.1.0)"]
 
 [[package]]
 name = "tomli"
-version = "2.0.1"
+version = "2.0.2"
 description = "A lil' TOML parser"
 optional = false
-python-versions = ">=3.7"
+python-versions = ">=3.8"
 files = [
-    {file = "tomli-2.0.1-py3-none-any.whl", hash = "sha256:939de3e7a6161af0c887ef91b7d41a53e7c5a1ca976325f429cb46ea9bc30ecc"},
-    {file = "tomli-2.0.1.tar.gz", hash = "sha256:de526c12914f0c550d15924c62d72abc48d6fe7364aa87328337a31007fe8a4f"},
+    {file = "tomli-2.0.2-py3-none-any.whl", hash = "sha256:2ebe24485c53d303f690b0ec092806a085f07af5a5aa1464f3931eec36caaa38"},
+    {file = "tomli-2.0.2.tar.gz", hash = "sha256:d46d457a85337051c36524bc5349dd91b1877838e2979ac5ced3e710ed8a60ed"},
 ]
 
 [[package]]
@@ -2100,4 +2100,4 @@ files = [
 [metadata]
 lock-version = "2.0"
 python-versions = ">=3.10,<4"
-content-hash = "1925239638206c6776ac6e50233e649e379ccb976212ecbea55be127d22665b2"
+content-hash = "c25e1b7093b0bbbc49c519efd3a87d5ccf63e7e78c899bf013e8181c9157cb3c"
diff --git a/pyproject.toml b/pyproject.toml
index 3dedf9cc9..9c4c883c0 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -32,7 +32,7 @@ pytest-cov = "5.0.0"
 coverage = "7.6.1"
 sympy = "1.13.3"
 distlib = "0.3.8"
-ipython = "8.27.0"
+ipython = "8.28.0"
 setuptools = ">65.5.1"
 
 [tool.poetry.group.lint.dependencies]
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_orthonormal.py b/toqito/matrix_props/is_orthonormal.py
index 876fb3695..244021853 100644
--- a/toqito/matrix_props/is_orthonormal.py
+++ b/toqito/matrix_props/is_orthonormal.py
@@ -39,4 +39,5 @@ def is_orthonormal(vectors: list[np.ndarray]) -> bool:
     :return: True if vectors are orthonormal; False otherwise.
 
     """
-    return is_mutually_orthogonal(vectors) and np.allclose(np.dot(vectors, vectors.T), np.eye(vectors.shape[0]))
+    return is_mutually_orthogonal(vectors) and np.allclose(np.dot(vectors, np.conjugate(vectors).T), np.eye(
+        vectors.shape[0]))
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/nonlocal_games/extended_nonlocal_game.py b/toqito/nonlocal_games/extended_nonlocal_game.py
index 933edf7e2..d778b62ac 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 7e42f0cfd..1f6dcc4f1 100644
--- a/toqito/nonlocal_games/nonlocal_game.py
+++ b/toqito/nonlocal_games/nonlocal_game.py
@@ -363,7 +363,6 @@ def __optimize_alice(self, dim, bob_povms) -> tuple[dict, float]:
         # .. math::
         #    \sum_{(x,y) \in \Sigma} \pi(x, y) V(a,b|x,y) \ip{B_b^y}{A_a^x}
         win = 0
-        is_real = True
         for x_ques in range(num_inputs_alice):
             for y_ques in range(num_inputs_bob):
                 for a_ans in range(num_outputs_alice):
@@ -372,17 +371,16 @@ 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],
                             cvxpy.expressions.variable.Variable,
                         ):
-                            is_real = False
                             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 +428,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/rand/__init__.py b/toqito/rand/__init__.py
index 9fd42e46d..94acf6fcd 100644
--- a/toqito/rand/__init__.py
+++ b/toqito/rand/__init__.py
@@ -7,3 +7,4 @@
 from toqito.rand.random_state_vector import random_state_vector
 from toqito.rand.random_states import random_states
 from toqito.rand.random_circulant_gram_matrix import random_circulant_gram_matrix
+from toqito.rand.random_orthonormal_basis import random_orthonormal_basis
diff --git a/toqito/rand/random_orthonormal_basis.py b/toqito/rand/random_orthonormal_basis.py
new file mode 100644
index 000000000..73df01d93
--- /dev/null
+++ b/toqito/rand/random_orthonormal_basis.py
@@ -0,0 +1,41 @@
+"""Generate random orthonormal basis."""
+import numpy as np
+
+from toqito.matrix_props import is_orthonormal
+from toqito.rand import random_unitary
+
+
+def random_orthonormal_basis(dim: int, is_real: bool = False) -> list[np.ndarray]:
+    r"""Generate a real random orthonormal basis of given dimension :math:`d`.
+
+    The basis is generated from the columns of a random unitary matrix of the same dimension
+    as the columns of a unitary matrix typically form an orthonormal basis :cite:`SE_1688950`.
+
+    Examples
+    ==========
+    To generate a random orthonormal basis of dimension :math:`4`,
+
+    >>> from toqito.rand import random_orthonormal_basis
+    >>> random_orthonormal_basis(4, is_real = True) # doctest: +SKIP
+    [array([0.52188745, 0.4983613 , 0.69049811, 0.04981832]),
+    array([-0.48670459,  0.58756912, -0.10226756,  0.63829658]),
+    array([ 0.23965404, -0.58538248,  0.187136  ,  0.75158061]),
+    array([ 0.658269  ,  0.25243989, -0.69118291,  0.158815  ])]
+
+    References
+    ==========
+    .. bibliography::
+        :filter: docname in docnames
+
+    dim: int
+        Number of elements in the random orthonormal basis.
+
+    """
+    random_mat = random_unitary(dim, is_real)
+
+    random_mat = random_unitary(dim, is_real)
+    return [random_mat[:, i] for i in range(dim)]
+
+    return rand_orth_basis
+
+


Test patch:
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/rand/tests/test_random_orthonormal_basis.py b/toqito/rand/tests/test_random_orthonormal_basis.py
new file mode 100644
index 000000000..cc9a0d419
--- /dev/null
+++ b/toqito/rand/tests/test_random_orthonormal_basis.py
@@ -0,0 +1,17 @@
+"""Tests for random orthonormal basis."""
+
+import numpy as np
+import pytest
+
+from toqito.matrix_props import is_orthonormal
+from toqito.rand import random_orthonormal_basis
+
+
+@pytest.mark.parametrize("input_dim", range(2, 5))
+@pytest.mark.parametrize("bool", [False, True])
+def test_random_orth_basis_int_dim(input_dim, bool):
+    """Test function works as expected for a valid int input."""
+    gen_basis = random_orthonormal_basis(dim=input_dim, is_real = bool)
+    assert len(gen_basis) == input_dim
+    assert is_orthonormal(np.array(gen_basis))
+