Skip to content

parameterize random_psd_operator with wishart option#1390

Merged
vprusso merged 22 commits intovprusso:masterfrom
aman-coder03:feature/random-psd-parameterization
Apr 17, 2026
Merged

parameterize random_psd_operator with wishart option#1390
vprusso merged 22 commits intovprusso:masterfrom
aman-coder03:feature/random-psd-parameterization

Conversation

@aman-coder03
Copy link
Copy Markdown
Contributor

@aman-coder03 aman-coder03 commented Jan 23, 2026

Description

closes #1275

this PR parameterizes the distribution used in random_psd_operator by introducing a new distribution argument. This allows users to select alternative sampling strategies when generating random positive semidefinite (PSD) operators
the default behavior remains unchanged to preserve backward compatibility

Changes

  • Added distribution: str = "uniform" parameter to random_psd_operator
  • Preserved existing behavior as the default ("uniform") for full backward compatibility
  • Added "wishart" option using Gaussian sampling (X @ X.conj().T) to generate PSD operators
  • Added validation for invalid distribution values
  • Added unit tests for Wishart sampling and invalid distribution handling

Checklist

  • Use ruff for errors related to code style and formatting
  • Verify all previous and newly added unit tests pass in pytest (targeted tests for modified module pass locally)
  • Check the documentation build does not lead to any failures

@codecov
Copy link
Copy Markdown

codecov Bot commented Jan 23, 2026

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 98.3%. Comparing base (b9865b8) to head (4386100).

Additional details and impacted files
@@          Coverage Diff           @@
##           master   #1390   +/-   ##
======================================
  Coverage    98.3%   98.3%           
======================================
  Files         202     202           
  Lines        5212    5221    +9     
  Branches     1196    1200    +4     
======================================
+ Hits         5124    5133    +9     
  Misses         46      46           
  Partials       42      42           

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@vprusso
Copy link
Copy Markdown
Owner

vprusso commented Jan 23, 2026

@aman-coder03 the code coverage went down with your changes. Please ensure you have tests that cover the areas that you modified. Code coverage report can be found here:

https://app.codecov.io/gh/vprusso/toqito/pull/1390

@aman-coder03
Copy link
Copy Markdown
Contributor Author

i have covered all the areas now @vprusso
thanks!

Comment on lines +82 to +111

def test_random_psd_operator_wishart():
"""Test Wishart distribution generates PSD matrix."""
mat = random_psd_operator(4, distribution="wishart")
assert mat.shape == (4, 4)
assert is_positive_semidefinite(mat)


def test_random_psd_operator_invalid_distribution():
"""Test invalid distribution raises ValueError."""
with pytest.raises(ValueError):
random_psd_operator(4, distribution="invalid")

def test_random_psd_operator_invalid_dim_type():
"""Test invalid dim type raises ValueError."""
with pytest.raises(ValueError):
random_psd_operator("4")


def test_random_psd_operator_invalid_dim_negative():
"""Test negative dim raises ValueError."""
with pytest.raises(ValueError):
random_psd_operator(-2)


def test_random_psd_operator_wishart_real_branch():
"""Test Wishart distribution with real sampling branch."""
mat = random_psd_operator(4, is_real=True, distribution="wishart")
assert mat.shape == (4, 4)
assert is_positive_semidefinite(mat)
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There should be a way to use the pytest.parameterize decorator to consolidate these tests so that they are not all written as separate functions. You can refer to how this is used elsewhere in the repo for context/reference.

Copy link
Copy Markdown
Contributor Author

@aman-coder03 aman-coder03 Jan 23, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

sure @vprusso i will check and update shortly!

@aman-coder03
Copy link
Copy Markdown
Contributor Author

@vprusso i have updated the tests, please have a look!

@pytest.mark.parametrize(
"dim, distribution",
[
("4", "uniform"), # invalid dim type
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please keep the comment style consistent with other tests, e.g.

Suggested change
("4", "uniform"), # invalid dim type
# Invalid dim type.
("4", "uniform"),

etc.

Comment thread toqito/rand/random_psd_operator.py Outdated
Comment on lines +113 to +115
raise ValueError(
"Invalid distribution. Supported options are 'uniform' and 'wishart'."
)
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Linter is set for 120 char-width:

Suggested change
raise ValueError(
"Invalid distribution. Supported options are 'uniform' and 'wishart'."
)
raise ValueError("Invalid distribution. Supported options are 'uniform' and 'wishart'.")

@aman-coder03
Copy link
Copy Markdown
Contributor Author

i have implemented the suggestions @vprusso

@Newtech66
Copy link
Copy Markdown
Contributor

Newtech66 commented Jan 24, 2026

Any reason for not allowing greater customizability in the "wishart" option? The Wishart distribution is parametrized by a scale matrix $V$ and a parameter $n$ as given here. The current implementation amounts to setting $V$ to the identity matrix and $n=$dim. Also, why not use scipy.linalg.wishart?

@vprusso
Copy link
Copy Markdown
Owner

vprusso commented Jan 24, 2026

Any reason for not allowing greater customizability in the "wishart" option? The Wishart distribution is parametrized by a scale matrix V and a parameter n as given here. The current implementation amounts to setting V to the identity matrix and n = dim. Also, why not use scipy.linalg.wishart?

@aman-coder03 FYI These are all good suggestions. Thanks, @Newtech66 !

@aman-coder03
Copy link
Copy Markdown
Contributor Author

aman-coder03 commented Jan 24, 2026

@Newtech66 @vprusso thanks for your suggestions!!
my initial goal was to introduce parameterization in a minimal and backward compatible way to address the eigenvalue spread concern, without expanding the function signature too much in a single PR
i agree that allowing users to specify the scale matrix and degrees of freedom would make the implementation more faithful to the general Wishart definition
Using scipy.stats.wishart is also a reasonable option
if there is interest, i would be happy to extend this further to support a fully parameterized Wishart distribution in a follow up PR

@aman-coder03
Copy link
Copy Markdown
Contributor Author

@vprusso if there are no additional changes, can we merge this PR?

@vprusso
Copy link
Copy Markdown
Owner

vprusso commented Feb 5, 2026

Any reason for not allowing greater customizability in the "wishart" option? The Wishart distribution is parametrized by a scale matrix V and a parameter n as given here. The current implementation amounts to setting V to the identity matrix and n = dim. Also, why not use scipy.linalg.wishart?

@aman-coder03 would it not be possible to address this comment as well, @aman-coder03 ?

@aman-coder03
Copy link
Copy Markdown
Contributor Author

You’re right that the current implementation corresponds to the identity-scale Wishart case with df = dim. My initial goal was to introduce distribution parameterization in a minimal and backward-compatible way.

That said, I agree that supporting optional scale and df parameters would make the implementation more complete. I’ll extend this PR to include a fully parameterized Wishart option with proper validation and tests.

Thanks again for the helpful feedback. I will push an update shortly

@aman-coder03 aman-coder03 requested a review from vprusso February 9, 2026 14:11
"is_real",
[True, False],
)
def test_random_psd_operator_wishart(is_real):
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not include this test to parameterize over the other distributions as well? Indeed, this could just be rolled into the initial test that's the first test in this file (instead of creating a whole new separate and fractured one here), right?

Comment thread toqito/rand/random_psd_operator.py Outdated
Comment on lines 69 to 71
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
seeded = random_psd_operator(2, is_real=True, seed=42)
seeded

Comment thread toqito/rand/random_psd_operator.py Outdated
Comment on lines 59 to 60
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
from toqito.matrix_props import is_positive_semidefinite
is_positive_semidefinite(real_psd_mat)

Comment thread toqito/rand/random_psd_operator.py Outdated
Comment on lines 50 to 52
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
real_psd_mat = random_psd_operator(2, is_real=True)
real_psd_mat

Comment thread toqito/rand/random_psd_operator.py Outdated
Comment on lines 30 to 32
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
complex_psd_mat = random_psd_operator(2)
complex_psd_mat

Comment thread toqito/rand/random_psd_operator.py Outdated
Comment on lines 16 to 18
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There isn't anything in this description about the uniform and wishart distributions.

Comment thread toqito/rand/random_psd_operator.py Outdated
if scale is None:
scale = np.eye(dim)

if df is None:
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

df is generally used to indicate a dataframe which is not what this is here.

Comment thread toqito/rand/random_psd_operator.py Outdated

rand_mat = (rand_mat.conj().T + rand_mat) / 2
eigenvals, eigenvecs = np.linalg.eigh(rand_mat)
Q, _ = np.linalg.qr(eigenvecs)
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Upper case letters are not appropriate for variable names according to PEP-8

@aman-coder03
Copy link
Copy Markdown
Contributor Author

hey @vprusso, addressed all the feedback! fixed the docstring format, added validations for dim, scale, and num_degrees, corrected the complex wishart scaling and added the missing tests. please have a look when you get a chance!


@pytest.mark.parametrize(
"dim, is_real, seed, expected_mat",
"dim, is_real, seed, distribution, expected_mat",
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is going on with the formatting here?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@vprusso that's a GitHub UI artifact from the merge conflict resolution, old and new parameter strings got rendered on the same line in the diff view. the actual code is correctly formatted. just verified locally with Select-String and all 20 tests pass, please have a look at the raw file!

Copy link
Copy Markdown
Owner

@vprusso vprusso left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the contribution! The core implementation looks solid and well-tested. A few things to address before we can merge:

Blocking

1. Docstring format regression

The PR uses Sphinx :code: and :math: directives, but the codebase has migrated to mkdocs/markdown. Compare with random_density_matrix.py which uses backticks for inline code (`dim`) and \(\) for inline math. The :code: and :math: syntax will render as raw text in the docs.

Fix: replace all :code:\"..."`with ``"..."`` and:math:\text{dim}=2with(\text{dim}=2)`.

2. Docstring examples: redundant imports and missing blank lines

The existing style (see random_density_matrix.py) puts blank lines between import, call, and print within each code block, and avoids re-importing modules already available in the same session. The Wishart example re-imports random_psd_operator despite being in session="psd_operator" where it's already imported.

Should Fix

3. Empty References section

    References:
        .. footbibliography::

This is a Sphinx directive that does nothing in mkdocs. Either add an actual citation to the Wishart distribution using the project's [@citation_key] / refs.bib system, or remove the section entirely.

4. num_degrees < dim produces singular matrices

When num_degrees < dim, the Wishart matrix is guaranteed to be rank-deficient. While technically valid PSD, this is likely unintentional for most users. At minimum, document this behavior in the docstring. A stronger option would be to emit a warning.

5. Wishart seeded tests should have hardcoded expected values

The uniform seeded tests assert against hardcoded expected matrices (good -- catches regressions). The Wishart seeded tests only check reproducibility (matrix == matrix2) but not against a fixed expected value. A regression that consistently produces wrong output would pass undetected. Please add hardcoded expected matrices for the Wishart seed cases too.

Nits

6. The non-PSD scale matrix test case is a ~140-char single line. Break it across lines for readability.

@aman-coder03
Copy link
Copy Markdown
Contributor Author

@vprusso i addressed all the feedback! fixed the docstring format (backticks and \(\) for math), removed the sphinx references section, added a warning for num_degrees < dim, hardcoded expected values for the wishart seeded tests, and cleaned up the non PSD scale matrix test. please have a look when you get a chance!

@aman-coder03 aman-coder03 requested a review from vprusso March 7, 2026 17:39
@purva-thakre purva-thakre force-pushed the master branch 4 times, most recently from 5fdc8f7 to 25f9e63 Compare March 18, 2026 07:26
Copy link
Copy Markdown
Owner

@vprusso vprusso left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for the work on this @aman-coder03. A few issues that need to be addressed before this can be merged:

Duplicate code fence markers

Lines 54, 62, and 74 each have a duplicate opening fence:

        ```python exec="1" source="above" session="psd_operator"
        ```python exec="1" source="above" result="text" session="psd_operator"

This will break the docs build. Only one opening fence per code block — remove the duplicate lines (the ones without result="text").

Args section is duplicated

The Args section appears twice in the docstring — once near the top (lines 27-34 in the original) and once in the new code (after line 81). The second one overrides the first, but the first one is now orphaned. Remove the first Args block and keep only the comprehensive one at the bottom.

dim validation is a behavior change

Adding if not isinstance(dim, int) or dim < 1: raise ValueError changes the existing API contract. Previously dim=0 or float dims would just produce empty/weird matrices. While the validation is reasonable, it should be called out as a breaking change. It also doesn't match other rand/ functions which don't validate dim this way (e.g., random_state_vector, random_density_matrix). I'd remove this validation for consistency — let numpy handle bad dims naturally.

Wishart complex scaling comment

The / np.sqrt(2) scaling comment says it ensures the expected scale matrix is scale rather than 2 * scale. This is correct for the standard complex Wishart definition, but worth adding a brief reference (e.g., "see Goodman (1963) complex Wishart distribution").

Test rtol change

The rtol=1e-6 change in assert_allclose for the seeded tests is fine, but should be explained — why was the default tolerance insufficient?

Minor

  • The Raises section in the docstring should list ValueError with descriptions of all the conditions that can trigger it (invalid dim, invalid distribution, bad scale, bad num_degrees).

@aman-coder03
Copy link
Copy Markdown
Contributor Author

@vprusso , addressed all the feedback, ruff checks clean and all 18 tests passing.

@aman-coder03 aman-coder03 requested a review from vprusso March 19, 2026 16:10
Copy link
Copy Markdown
Owner

@vprusso vprusso left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for picking up #1275 — the overall shape of this is good: backward compatibility is preserved (existing seeded uniform expected matrices are unchanged, positional call random_psd_operator(dim, is_real, seed) still works), the docstring got a nice cleanup, and the Wishart normalization convention with the /sqrt(2) factor and Goodman (1963) citation is correct. A handful of things to address before merging.

Should fix

  1. scale must survive list/tuple input or the type annotation should forbid it. scale: np.ndarray | None is the declared type, but most users will pass a plain list. The current code does scale.shape != (dim, dim) directly, which crashes with AttributeError instead of a clean ValueError. Add scale = np.asarray(scale) as the first line of the else branch (or coerce on entry).

  2. num_degrees is not checked for being an integer. A caller passing num_degrees=2.5 passes the < 1 check and then crashes deep inside np.random.Generator.multivariate_normal with an unhelpful error. Add isinstance(num_degrees, (int, np.integer)) to the validation.

  3. test_random_psd_operator_invalid_distribution doesn't use match=. Right now any ValueError anywhere in the function satisfies this test — e.g., if a future refactor makes it raise for a different reason, the test still passes. Please match on "Supported options" or similar so the test actually pins the error path.

Style / clarity

  1. Complex Wishart draws multivariate_normal twice. This works but is wasteful and a little awkward. A cleaner alternative is to draw once at double width and split, or (standard) Cholesky-factor the scale matrix and multiply a white complex Gaussian:

    L = np.linalg.cholesky(scale)  # fall back to eigendecomp if scale is singular
    z = (gen.standard_normal((dim, num_degrees))
         + 1j * gen.standard_normal((dim, num_degrees))) / np.sqrt(2)
    x_mat = L @ z

    Same distribution, half the MVN calls, and it composes naturally if scale is later allowed to be rank-deficient. Not blocking — current code is correct.

  2. distribution could be Literal["uniform", "wishart"] for better type-checker hints and autocomplete, rather than bare str.

  3. rtol=1e-6 comment says "due to floating-point accumulation in Wishart matrix construction" — if the same seed + same code should produce bit-identical output on the same machine, the likely real reason you needed a looser tolerance is cross-platform BLAS nondeterminism. Consider updating the comment to say that, otherwise it reads like the math is intrinsically lossy.

  4. scale is not None or num_degrees is not None on the uniform path emits a warning. Reasonable, but arguably this should be a ValueError — silently ignoring parameters the user explicitly passed is a footgun, and since there are only two distributions right now the user's intent is unambiguous. Judgment call; happy either way but worth a thought.

Out of scope (just flagging)

  1. The pre-existing uniform path does q_mat, _ = np.linalg.qr(eigenvecs) on already-unitary eigenvectors from np.linalg.eigh, which is redundant — the QR of a unitary matrix is itself up to a diagonal sign. This isn't introduced by your PR, so don't fix it here, but worth a follow-up issue to simplify the uniform branch to eigenvecs @ diag(|eigvals|) @ eigenvecs.conj().T.

if scale is None:
scale = np.eye(dim)
else:
if scale.shape != (dim, dim):
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

A caller passing a plain list (e.g., scale=[[1,0],[0,1]]) will crash here with AttributeError: 'list' object has no attribute 'shape' instead of the clean ValueError you're aiming for — and the type annotation scale: np.ndarray | None doesn't actually prevent that at runtime.

Suggest coercing first:

else:
    scale = np.asarray(scale)
    if scale.shape != (dim, dim):
        raise ValueError(f"scale must be a {dim}x{dim} matrix, got {scale.shape}.")
    if not is_positive_semidefinite(scale):
        raise ValueError("scale must be a positive semidefinite matrix.")

Comment thread toqito/rand/random_psd_operator.py Outdated

if num_degrees is None:
num_degrees = dim
if num_degrees < 1:
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The error message says "must be a positive integer" but you only check the positivity, not the integer-ness. num_degrees=2.5 slips through here and then crashes inside np.random.Generator.multivariate_normal with a less helpful error.

if not isinstance(num_degrees, (int, np.integer)) or num_degrees < 1:
    raise ValueError("num_degrees must be a positive integer.")

# based on a certain multivariate complex Gaussian distribution."
x_mat = (
gen.multivariate_normal(np.zeros(dim), scale, size=num_degrees).T
+ 1j * gen.multivariate_normal(np.zeros(dim), scale, size=num_degrees).T
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Minor: this calls multivariate_normal twice per invocation. You can generate the complex Wishart sample more directly by drawing a single white complex Gaussian and left-multiplying by the Cholesky factor of scale:

L = np.linalg.cholesky(scale)
z = (gen.standard_normal((dim, num_degrees))
     + 1j * gen.standard_normal((dim, num_degrees))) / np.sqrt(2)
x_mat = L @ z

Same distribution (E[x_mat x_mat^H] = num_degrees * scale), one RNG path, and it keeps the real and complex branches structurally similar. Not blocking — current code is correct.

) / np.sqrt(2)
return x_mat @ x_mat.conj().T

raise ValueError("Invalid distribution. Supported options are 'uniform' and 'wishart'.")
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Consider typing distribution as Literal["uniform", "wishart"] in the signature so type checkers and IDE autocomplete can surface the valid options. The runtime ValueError here is still the right safety net for untyped callers.

assert_allclose(matrix, expected_mat, rtol=1e-6)


def test_random_psd_operator_invalid_distribution():
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

pytest.raises(ValueError) without match= is too permissive — any ValueError raised anywhere in the function will satisfy this assertion, including ones raised for unrelated reasons if the code is later refactored. Please add match="Supported options" (or similar) so the test actually pins the intended error path.

with pytest.raises(ValueError, match="Supported options"):
    random_psd_operator(4, distribution="invalid")

assert_allclose(matrix, expected_mat)
# rtol=1e-6 used due to floating-point accumulation in Wishart matrix construction.
matrix = random_psd_operator(dim, is_real, seed, distribution=distribution)
assert_allclose(matrix, expected_mat, rtol=1e-6)
Copy link
Copy Markdown
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Nit: the inline comment attributes rtol=1e-6 to "floating-point accumulation in Wishart matrix construction," but if the same seed + same NumPy version produce the same output on your machine to machine precision, the real reason you need a looser tolerance is more likely cross-platform BLAS/LAPACK nondeterminism. Worth updating the comment so future maintainers understand why the tolerance is loose.

@aman-coder03 aman-coder03 requested a review from vprusso April 11, 2026 14:01
@aman-coder03
Copy link
Copy Markdown
Contributor Author

you can have a look now @vprusso

@vprusso vprusso merged commit 2e6c749 into vprusso:master Apr 17, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Eigenvalue spread of random_psd_operator is very skewed

3 participants