Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Matrix creation from a scalar fails in some cases #36065

Closed
2 tasks done
vneiger opened this issue Aug 10, 2023 · 1 comment
Closed
2 tasks done

Matrix creation from a scalar fails in some cases #36065

vneiger opened this issue Aug 10, 2023 · 1 comment
Labels

Comments

@vneiger
Copy link
Contributor

vneiger commented Aug 10, 2023

Steps To Reproduce

Creating a matrix from a scalar (MA_ENTRIES_SCALAR flag from args.pyx), using code such as mat = matrix(ZZ, nrows, ncols, value), will unexpectedly fail in some cases. Such a call is supposed to build a diagonal matrix with value on the diagonal, and requires the matrix to be square (nrows == ncols) if value is nonzero. The issue is mainly that we cannot always successfully compare value to zero...

Consider the following example, inspired from a test in the doc for get_is_zero_unsafe in matrix0.pyx :

sage: class MyAlgebraicNumber(sage.rings.qqbar.AlgebraicNumber):
....:     def __bool__(self):
....:         raise ValueError
....:
sage: mat = matrix(1,1,MyAlgebraicNumber(1))
sage: mat
[1]
sage: mat = matrix(1,2,MyAlgebraicNumber(1))
### FAILS with ValueError ###

The failure is the same if using MyAlgebraicNumber(0).

Expected Behavior

The wanted behaviour is that such a call fails when nrows != ncols and value is not zero. The current error, for example if one does matrix(ZZ, 1, 2, 1) is TypeError("nonzero scalar matrix must be square"). If value happens to be zero, then the code just creates the zero matrix (even when nrows != ncols).

Actual Behavior

We get a value error, as below.

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In [4], line 1
----> 1 mat = matrix(Integer(1),Integer(2),MyAlgebraicNumber(Integer(1)))

File ~/repositories/software/sage/src/sage/matrix/constructor.pyx:643, in sage.matrix.constructor.matrix()
    641 """
    642 immutable = kwds.pop('immutable', False)
--> 643 M = MatrixArgs(*args, **kwds).matrix()
    644 if immutable:
    645     M.set_immutable()

File ~/repositories/software/sage/src/sage/matrix/args.pyx:655, in sage.matrix.args.MatrixArgs.matrix()
    653     True
    654 """
--> 655 self.finalize()
    656
    657 cdef Matrix M

File ~/repositories/software/sage/src/sage/matrix/args.pyx:929, in sage.matrix.args.MatrixArgs.finalize()
    927 if self.typ == MA_ENTRIES_SCALAR:
    928     if self.nrows != self.ncols:
--> 929         if self.entries:
    930             raise TypeError("nonzero scalar matrix must be square")
    931         self.typ = MA_ENTRIES_ZERO

Cell In [1], line 3, in MyAlgebraicNumber.__bool__(self)
      2 def __bool__(self):
----> 3     raise ValueError

ValueError:

Additional Information

The culprit is the piece of finalize in args.pyx as highlighted in the above error message. I would suggest to correct the code so that it fails if the matrix is rectangular and value cannot be determined to be zero. So, on the above examples with both and , it would fail with a TypeError with a message such as "a scalar matrix must be square if it is nonzero or cannot be determined to be zero".

Environment

- **OS**: Ubuntu 22.04
- **Sage Version**: 10.1.beta9

Checklist

  • I have searched the existing issues for a bug report that matches the one I want to file, without success.
  • I have read the documentation and troubleshoot guide
@vneiger vneiger added the t: bug label Aug 10, 2023
vneiger added a commit to vneiger/sage that referenced this issue Aug 10, 2023
vbraun pushed a commit to vbraun/sage that referenced this issue Aug 12, 2023
… type MA_ENTRIES_ZERO

    
There are many ways to specify entries when creating a matrix, which are
handled in `args.pyx` with `MatrixArgs` class. It has been discussed
before why allowing `entries=None` in the matrix creation methods can be
beneficial in terms of performance
(sagemath#11589 ,
sagemath#12020). This input
`entries=None` is required to yield the zero matrix. For example:
`matrix(QQ, 4, 4)` creates the zero matrix. This corresponds to the type
`MA_ENTRIES_ZERO` in `args.pyx`. When one passes a value, e.g.
`matrix(QQ, 4, 4, 2)`, then one gets a diagonal matrix with `2` on the
diagonal; this corresponds to the type `MA_ENTRIES_SCALAR` in
`args.pyx`.

Currently, doing something like `matrix(QQ, 4, 4, 0)` will pick
`MA_ENTRIES_SCALAR`, and therefore will build the matrix and fill the
diagonal with zero. [Behind the scenes, there is still some
acknowledgement that this is not the usual scalar matrix case, since
this will not fail if the matrix is not square (`matrix(QQ, 4, 5, 0)`
will not fail, but `matrix(QQ, 4, 5, 1)` will). But this is still not
seen as `MA_ENTRIES_ZERO`.] This PR ensures the type `MA_ENTRIES_ZERO`
is picked in this situation. This improves performance and solves some
issues, as noted below.  This PR also fixes the related issue
sagemath#36065 .

In fact, `entries=None` is the default in the `__init__` of all matrix
subclasses presently implemented. It also used to be the default when
constructing a matrix by "calling" a matrix space, until
sagemath#31078 and https://github.com/sag
emath/sage/commit/cef613a0a57b85c1ebc5747185213ae4f5ec35f2 which changed
this default value from `None` to `0`, bringing some performance
regression.

Regarding this "calling the matrix space", this PR also solves the
performance issue noted in
sagemath#35961 (comment) ,
where it was observed that in the following piece of code:
```
sage: space = MatrixSpace(GF(9001), 10000, 10000)
sage: %timeit zmat = space.matrix()
18.3 µs ± 130 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops
each)
sage: %timeit zmat = space()
12.1 ms ± 65.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
```
the called default is not the same. `space.matrix()` directly
initializes the matrix through `entries=None`, but on the other hand
`space()` actually calls the constructor with `entries=0`. This
performance regression comes from
sagemath#31078 and https://github.com/sag
emath/sage/commit/cef613a0a57b85c1ebc5747185213ae4f5ec35f2 , where the
default for construction from the matrix space was changed from `None`
to `0`. This cannot be easily reverted: this is now the default in the
`__call__` of the `Parent` class. So this PR does not directly revert
the call default to `None`  somehow, but the result is very close in
effect (read: the overhead is small). Unchanged: through the parent
call, `0` is converted to the base ring zero, which is passed to the
constructor's `entries` which is then analyzed as `MA_ENTRIES_SCALAR`.
Changed: the code modifications ensure that soon enough it will be
detected that this is in fact `MA_ENTRIES_ZERO`. The small overhead
explains why, despite the improvements, construction with `None` is
still sometimes slightly faster than with `0`.

Below are some timings showing the improvements for some fields. Also,
this PR merged with the improvements discussed in
sagemath#35961 will make the above
timings of `space.matrix()` and `space()` be the same (which means a
speed-up of a factor >500 for this call `space()`...). The measurement
is for construction via calling the matrix space: `None` is
`space(None)`, `Empty` is `space()`, `Zero` is `space(0)`, and `Nonzero`
is `space(some nonzero value)`.

```
NEW TIMES

field   dim     None    Empty   Zero    Nonzero
GF(2)   5       2.3e-06 3.2e-06 3.6e-06 3.2e-06
GF(2)   50      2.4e-06 3.3e-06 3.6e-06 5.8e-06
GF(2)   500     3.6e-06 4.5e-06 4.8e-06 3.1e-05
GF(512) 5       2.6e-05 2.8e-05 2.9e-05 2.9e-05
GF(512) 50      2.6e-05 2.9e-05 2.9e-05 4.0e-05
GF(512) 500     3.7e-05 3.8e-05 3.9e-05 1.6e-04
QQ      5       2.2e-06 3.3e-06 3.4e-06 3.2e-06
QQ      50      8.0e-06 9.2e-06 9.4e-06 1.2e-05
QQ      500     6.1e-04 6.3e-04 6.4e-04 6.7e-04

OLD TIMES

field   dim     None    Empty   Zero    Nonzero
GF(2)   5       2.3e-06 3.5e-06 3.6e-06 3.7e-06
GF(2)   50      2.4e-06 6.0e-06 6.1e-06 6.0e-06
GF(2)   500     3.6e-06 3.0e-05 3.0e-05 3.0e-05
GF(512) 5       2.5e-05 2.8e-05 2.9e-05 2.9e-05
GF(512) 50      2.5e-05 3.9e-05 4.0e-05 4.0e-05
GF(512) 500     3.5e-05 1.5e-04 1.5e-04 1.6e-04
QQ      5       2.2e-06 3.5e-06 3.7e-06 3.7e-06
QQ      50      7.9e-06 1.2e-05 1.2e-05 1.2e-05
QQ      500     6.4e-04 6.9e-04 6.9e-04 6.9e-04
```

Code used for the timings:
```
time_kwds = dict(seconds=True, number=20000, repeat=7)

fields = [GF(2), GF(2**9), QQ]
names = ["GF(2)", "GF(512)", "QQ"]
vals = [GF(2)(1), GF(2**9).gen(), 5/2]

print(f"field\tdim\tNone\tEmpty\tZero\tNonzero")
for field,name,val in zip(fields,names,vals):
    for dim in [5, 50, 500]:
        space = MatrixSpace(field, dim, dim)
        tnone = timeit("mat = space(None)", **time_kwds)
        tempty = timeit("mat = space()", **time_kwds)
        tzero = timeit("mat = space(0)", **time_kwds)
        tnonz = timeit("mat = space(1)", **time_kwds)
        print(f"{name}\t{dim}\t{tnone:.1e}\t{tempty:.1e}\t{tzero:.1e}\t{
tnonz:.1e}")
```

### 📝 Checklist

- [x] The title is concise, informative, and self-explanatory.
- [x] The description explains in detail what this PR is about.
- [x] I have linked a relevant issue or discussion.
- [x] I have created tests covering the changes.
- [x] I have updated the documentation accordingly.

### ⌛ Dependencies
    
URL: sagemath#36068
Reported by: Vincent Neiger
Reviewer(s): Matthias Köppe, Vincent Neiger
@vneiger
Copy link
Contributor Author

vneiger commented Aug 16, 2023

Fixed via #36068

@vneiger vneiger closed this as completed Aug 16, 2023
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
Projects
None yet
Development

No branches or pull requests

1 participant