random_matrix() with prescribed density buggy

[sagetrac-rpw]

Matrices with prescribed density are not generated correctly:

sage: M = random_matrix(GF(65537), 100, 100, sparse=True, density=0.1)
sage: len(M.nonzero_positions())
940
sage: M = random_matrix(GF(2), 100, 100, sparse=True, density=0.1)
sage: len(M.nonzero_positions())
465

To wit: the actual density of the matrix over GF(2) is only approximately half of what we expect. This happens because the randomize() function populating the entries does not check whether the random element picked actually is non-zero. Apparently, all of the matrix classes are affected by this bug.

CC: @craigcitro

Component: linear algebra

Author: Sebastian Pancratz

Reviewer: Tom Boothby, Craig Citro

Merged: sage-4.3.2.alpha0

Issue created by migration from https://trac.sagemath.org/ticket/3436

[jhpalmieri]

comment:2

I'll also point out that the documentation for random_matrix says that density should be an integer, not a float -- see also randomize in matrix2.pyx and random_element in matrix_space.py -- while in matrix_integer_dense.pyx, randomize says that density should be a float.

[jasongrout]

comment:3

Some more examples with sage 4.1. The last few examples show that the density is about an order of magnitude smaller than what we asked for.

sage: A=random_matrix(ZZ,1000,density=1e-3,sparse=True)  
sage: RR(len(A.nonzero_positions()))/(A.nrows()*A.ncols())
0.000814000000000000
sage: A=random_matrix(ZZ,10000,density=1e-4,sparse=True)  
sage: RR(len(A.nonzero_positions()))/(A.nrows()*A.ncols())
0.0000796700000000000
sage: A=random_matrix(ZZ,100000,density=1e-5,sparse=True) 
sage: RR(len(A.nonzero_positions()))/(A.nrows()*A.ncols())
1.13230000000000e-6
sage: A=random_matrix(ZZ,100000,density=1e-6,sparse=True) 
sage: RR(len(A.nonzero_positions()))/(A.nrows()*A.ncols())
1.12600000000000e-7
sage: A=random_matrix(ZZ,100000,density=1e-7,sparse=True) 
sage: RR(len(A.nonzero_positions()))/(A.nrows()*A.ncols())
1.08000000000000e-8

[jhpalmieri]

comment:4

Replying to @jasongrout:

Some more examples with sage 4.1. The last few examples show that the density is about an order of magnitude smaller than what we asked for.

When I tried this, I got answers like the first two: the ratio of (fraction of nonzero entries) / (density) was about 4/5. This is consistent with the observation in the ticket description that the randomize function doesn't check whether the random element is zero: for integers, the documentation for random_element says that the default distribution gives Pr(X = 0) = 1/5.

I would try to write a patch for this but pretty much all of the matrix code is written in Cython, which I don't know, so anything I write would likely be buggy and slow things down by a factor of 100.

[sagetrac-spancratz]

comment:5

Attached a patch. It introduces a new method in ring.pyx, _random_nonzero_element, which by default calls random_element until a non-zero element has been obtained.

On the matrix end, I've introduced a new parameter nonzero to the randomize methods in the various matrix implementations. If this is set to True, the modified entries are guaranteed to be non-zero, whereas if it is set to False they are not. The method random_element on a matrix space now calls randomize on a zero matrix with nonzero=True.

Finally, this is now also used to ensure that generating a fraction field element does not cause a division by zero error.

Sebastian

[sagetrac-spancratz]

Author: spancratz

[sagetrac-spancratz]

Attachment: trac3436_a.patch.gz

First patch

[sagetrac-spancratz]

comment:6

I've attached the patch now. I've called it "First patch" since I've only doctested sage/matrix/ and adjusted some docstrings as needed, but I expect that a few more such little changes might come up during a full test.

Sebastian

[sagetrac-spancratz]

Attachment: trac3436_b.patch.gz

Additional patch (more doctests)

[sagetrac-spancratz]

comment:7

The additional patch takes care of the remaining three doctest failures. This should now pass all doctests done with ./sage -t devel/sage/sage.

Sebastian