Add random_bistochastic_matrices() special constructor

[JosePisco]

I add the matrix method random_bistochastic_matrix() that generates a random bistochastic matrix in the matrix space in input, provided that it is a subfield of the real numbers.

The method used squares each element of a random unitary matrix since it will necessarily give a bistochastic matrix.

You can call the function in two ways:

sage: from sage.matrix.constructor import random_bistochastic_matrix
sage: MS = MatrixSpace(QQ, 4) 
sage: B = random_bistochastic_matrix(MS)

or with

sage: MS = MatrixSpace(QQ, 4) 
sage: B = matrix.random_bistochastic(MS)

Note that over the real numbers, we have inexact matrices and .is_bistochastic() may fail.

sage: MS = MatrixSpace(RR, 4) 
sage: B = matrix.random_bistochastic(MS)
sage: B.is_bistochastic()  # random
False

📝 Checklist

• The title is concise and informative.
• The description explains in detail what this PR is about.
• I have linked a relevant issue or discussion.
• I have created tests covering the changes.
• I have updated the documentation and checked the documentation preview.

[github-actions]

Documentation preview for this PR (built with commit cb089f0; changes) is ready! 🎉
This preview will update shortly after each push to this PR.

[orlitzky]

You'll need a double-colon in TESTS:: so that these are considered a code block. You may also want to separate the independent examples with ::, like

TESTS::

    sage: example1

::

    sage: example2

It's academic since the tests aren't usually rendered to the documentation, but each :: creates a separate code block.

[orlitzky]

The is_subring(RR) approach will miss some cases unfortunately. The best "is real" check I've found is the one in random_unitary_matrix. Something like,

from sage.rings.real_lazy import RLF
F = self.base_ring()
if not (RLF.has_coerce_map_from(F) or F.has_coerce_map_from(RLF)):
    # does not look like a subring of the reals

[JosePisco]

What kind of cases do we miss by doing is_subring(RR) ?

[orlitzky]

I don't remember all of them, but the real ball/interval fields are one example:

sage: RBF
Real ball field with 53 bits of precision
sage: RBF.is_subring(RR)
False
sage: RBF.has_coerce_map_from(RLF)
True

sage: RIF
Real Interval Field with 53 bits of precision
sage: RIF.is_subring(RR)
False
sage: RIF.has_coerce_map_from(RLF)
True

Neither one is perfect though.

[JosePisco]

I'm not sure I understand all the details of sage's different real number fields and their coercion map but that sounds reasonable given that's what is being used in other matrix methods.

[orlitzky]

I think we can get every bistochastic matrix using Birkhoff's theorem. The (finite) set of permutation matrices are readily available, so all you'd need to do is generate a random convex combination, i.e. a list of positive coefficients that sum to one.

[JosePisco]

We could but that would require factorial complexity. Or we can sample $(n-1)^2 + 1$ coefficients, as in Birkhoff's algorithm, but then the method to construct the matrix (programming-wise) would require more code because we would also have to sample for random permutations and so on and this would be much less sustainable to maintain in the future.

The unitary approach is very reasonnable and I don't think missing on some bistochastic matrices is a big deal.

[orlitzky]

Taking a random sample of the permutations is actually a one-liner:

sage: perms = sorted(Permutations(5))
sage: from random import sample
sage: sample(perms, 10)
[[2, 1, 3, 5, 4],
 [4, 5, 2, 1, 3],
 [3, 4, 1, 5, 2],
 [5, 2, 4, 1, 3],
 [3, 1, 4, 2, 5],
 [3, 2, 1, 5, 4],
 [2, 4, 3, 5, 1],
 [3, 1, 4, 5, 2],
 [3, 1, 2, 4, 5],
 [5, 3, 1, 2, 4]]

I think it would be cool to compare the result using bistochastic_as_sum_of_permutations, but I won't insist on this. You took the time to write it, so you're the boss.

[orlitzky]

You actually don't want the double-colon here, because the next line ("Matrices with rational entries") is not part of the code block.

[JosePisco]

Ah got it, thanks! Sorry for the poor mistakes, I really appreciate you pointing them out. I'm all for strictly following style and rules.