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

Improve Matrix_mod2_dense solve_right performance #38843

Open
wants to merge 7 commits into
base: develop
Choose a base branch
from

Conversation

maple3142
Copy link

Use M4RI builtin mzd_solve_left function to solve equation.

Currently, using solve_right on large GF(2) matrix is pretty slow because it use a generic solving function in matrix2.pyx, but M4RI provides a mzd_solve_left function that does the same thing.

For example, try running this in sage's IPython REPL:

set_random_seed(12345)
n, m = 20000, 19968
A = random_matrix(GF(2), n, m)
x = random_vector(GF(2), m)
B = A*x
%time sol = A.solve_right(B)
assert A * sol == B

without this PR, it took

CPU times: user 1min 2s, sys: 190 ms, total: 1min 2s
Wall time: 15.2 s

with this PR, it becomes

CPU times: user 12.5 s, sys: 40.8 ms, total: 12.6 s
Wall time: 1.09 s

📝 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.

⌛ Dependencies

Use M4RI builtin `mzd_solve_left` function to solve equation.
Copy link

github-actions bot commented Oct 23, 2024

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

# although it is called mzd_solve_left, it does the same thing as solve_right
ret = mzd_solve_left(lhs, rhs, 0, check)
sig_off()
mzd_free(lhs)
Copy link
Contributor

Choose a reason for hiding this comment

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

Wouldn't this leak the lhs and rhs if user press ctrl+C?

Copy link
Author

Choose a reason for hiding this comment

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

You are right, should I wrap it with a try...finally block or just disable the ability to Ctrl-C?

Copy link
Author

Choose a reason for hiding this comment

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

Similarity, while _right_kernel_matrix does not enable signal but it can still leak the resulting kernel matrix when Ctrl-C (triggered during self.new_matrix).

Copy link
Contributor

Choose a reason for hiding this comment

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

Certainly handling it properly is preferable. (But no memory leak is a hard thing to test for…)

cysignals is a bit complicated…

https://cysignals.readthedocs.io/en/latest/interrupt.html#using-sig-on-and-sig-off

So far I don't exactly understand what's the difference between try: sig_on() except KeyboardInterrupt: versus if not sig_on_no_except(): either.

Copy link
Contributor

Choose a reason for hiding this comment

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

for what it's worth, the part I don't understand is sagemath/cysignals#205 .

In any case, I think the following will work:

data = malloc(...)
if not sig_on_no_except():
	mzd_free(lhs)
	mzd_free(rhs)
	cython_check_exception()

ret = mzd_solve_left(...)
sig_off()

Copy link
Author

Choose a reason for hiding this comment

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

Okay, I have pushed a version where try...finally are used to ensure allocated matrices get freed.

if self._entries.ncols == 0 or B_entries.ncols == 0:
# special case: empty matrix
if B != 0:
raise ValueError("matrix equation has no solutions")
Copy link
Contributor

Choose a reason for hiding this comment

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

Surely this is redundant (?) if B has no entries then every entries are 0?

(also even if it can happen you forget to free X here)

Copy link
Author

@maple3142 maple3142 Oct 24, 2024

Choose a reason for hiding this comment

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

There are cases where A have zero cols but B have non zero cols, so it would have a different result depends on whether B is zero. And without this, some sage tests will fail.

As for X, I think it is a Python object so it would be handled by Python's refcount?

@user202729
Copy link
Contributor

user202729 commented Oct 28, 2024

Maybe add a few tests?

sage: matrix(GF(2), 2, 0).solve_right(matrix(GF(2), 2, 2))
[]
sage: matrix(GF(2), 2, 0).solve_right(matrix(GF(2), 2, 2)+1)
Traceback…

ValueError: matrix equation has no solutions
sage: matrix(GF(2), 2, 0).solve_right(matrix(GF(2), 2, 2)+1, check=False)
[]
sage: matrix(GF(2), 2, 2).solve_right(matrix(GF(2), 2, 0))
[]

(note that the old code returns empty matrix for the third case, while the new code raises error. Apparently it doesn't break any existing test case, but anyway.)

Strange that m4ri doesn't automatically handle the empty matrix case correctly (mathematically) though. Any idea why?

@maple3142
Copy link
Author

Tests for special case added.

Strange that m4ri doesn't automatically handle the empty matrix case correctly (mathematically) though. Any idea why?

TBH, zero size matrix is pretty weird by itself I think, so it is not unusual to not consider that I think. It would probably be better handled by the solve_right in matrix2.py so the underlying function don't need to handle that.

@maple3142
Copy link
Author

maple3142 commented Nov 5, 2024

Not sure what breaks the test using Meson on Python 3.9...
The failing test does not fail locally for me too.

@user202729
Copy link
Contributor

user202729 commented Nov 5, 2024

My argument: just as zero is not a weird integer, zero-dimensional vector space (which is just a point) is not a weird space. So a zero-sized matrix is just a linear map from (or to) a zero-dimensional vector space.

The math checks out.

(Of course, programming-wise it is still a corner case.)


The failed test case (commit 9c532c5) is

2024-11-04T17:56:52.3683952Z File "src/sage/categories/simplicial_sets.py", line 460, in sage.categories.simplicial_sets.SimplicialSets.Pointed.ParentMethods.covering_map
2024-11-04T17:56:52.3685700Z Failed example:
2024-11-04T17:56:52.3686384Z     C.domain().face_data()
2024-11-04T17:56:52.3686928Z Expected:
2024-11-04T17:56:52.3687318Z     {(*, ()): None,
2024-11-04T17:56:52.3687752Z      (*, (1,2,3)): None,
2024-11-04T17:56:52.3688211Z      (*, (1,3,2)): None,
2024-11-04T17:56:52.3750815Z      (sigma_1, ()): ((*, (1,2,3)), (*, ())),
2024-11-04T17:56:52.3751444Z      (sigma_1, ()): ((*, ()), (*, ())),
2024-11-04T17:56:52.3752076Z      (sigma_1, (1,2,3)): ((*, (1,3,2)), (*, (1,2,3))),
2024-11-04T17:56:52.3752801Z      (sigma_1, (1,2,3)): ((*, (1,2,3)), (*, (1,2,3))),
2024-11-04T17:56:52.3753508Z      (sigma_1, (1,3,2)): ((*, ()), (*, (1,3,2))),
2024-11-04T17:56:52.3754218Z      (sigma_1, (1,3,2)): ((*, (1,3,2)), (*, (1,3,2)))}
2024-11-04T17:56:52.3754795Z Got:
2024-11-04T17:56:52.3755155Z     {(*, ()): None,
2024-11-04T17:56:52.3755553Z      (*, (1,2,3)): None,
2024-11-04T17:56:52.3755990Z      (*, (1,3,2)): None,
2024-11-04T17:56:52.3756430Z      (sigma_1, ()): ((*, ()), (*, ())),
2024-11-04T17:56:52.3757008Z      (sigma_1, ()): ((*, (1,2,3)), (*, ())),
2024-11-04T17:56:52.3757684Z      (sigma_1, (1,2,3)): ((*, (1,2,3)), (*, (1,2,3))),
2024-11-04T17:56:52.3758355Z      (sigma_1, (1,2,3)): ((*, (1,3,2)), (*, (1,2,3))),
2024-11-04T17:56:52.3759027Z      (sigma_1, (1,3,2)): ((*, (1,3,2)), (*, (1,3,2))),
2024-11-04T17:56:52.3759724Z      (sigma_1, (1,3,2)): ((*, ()), (*, (1,3,2)))}

looks completely unrelated to this change. And the only difference is in the ordering of the keys so… probably it's fine.

On the other hand Python dict is ordered by insertion time, and it looks weird that a dict can have two keys that is exactly equal…? (they're probably not equal but have same representation or something. Worth opening a new issue for that.)

@user202729
Copy link
Contributor

Looks like the failing test is #38940 . Should be unrelated to this change then.

sage: matrix(GF(2), 2, 0).solve_right(matrix(GF(2), 2, 2) + 1, check=False) == matrix(GF(2), 0, 2)
True
sage: matrix(GF(2), 2, 2).solve_right(matrix(GF(2), 2, 0)) == matrix(GF(2), 2, 0)
True
Copy link
Contributor

Choose a reason for hiding this comment

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

Suggested change
True
True
Check that it can be interrupted::
sage: set_random_seed(12345)
sage: n, m = 20000, 19968
sage: A = random_matrix(GF(2), n, m)
sage: x = random_vector(GF(2), m)
sage: B = A*x
sage: alarm(0.5); sol = A.solve_right(B)
Traceback (most recent call last):
...
AlarmInterrupt

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

2 participants