# Linear Algebra over different fields (`LinearAlgebra.det` does not respect the `+` and `*` multimethods)

I want to get a random matrix of binary numbers (the GF(2) field, i.e. `{0,1}`), then check whether it is invertible by calculating its determinant and getting its inverse.

From stylistic (and elegance of types) point of view, I was hoping that the following would work:

``````using LinearAlgebra
bin_matrix = rand(Bool, 4, 4)
determinant = det(bin_matrix)
if determinant != false
inverted = inv(bin_matrix)
end
``````

I expected `determinant` to be `Bool`, but it was `Float`. I also expected `inverted` to contain only `{true, false}`, but it actually contained `{0.0, 1.0, -1.0}`.

I know the algorithm for the correct determinant and inversion functions in principle, this is not my question.

Rather, I was wondering from the point of view of â€śJuliaâ€™s philosophyâ€ť, why are `Bool` matrices forced to behave as `Float` matrices?

It is especially surprising for me, because `true+false` and `true*false` in Julia do behave exactly like I would have expected (and I will be using them in my implementation of `det`), but it seems they are not the operations actually used in the computation of the `LinearAlgebra.det`? Isnâ€™t the whole point of the elegant multiple dispatch paradigm that functions like `LinearAlgebra.det` would actually use the special `+` and `*` methods that `Bool` provides?

Juliaâ€™s `LinearAlgebra` library doesnâ€™t have built-in functions for exact determinants of integer matrices, which is more appropriate for a number-theory package. Instead, its `det` function works by first computing the LU factorization (i.e. doing Gaussian elimination), which requires conversion to floating point (and is accurate up to roundoff errors).

3 Likes

Thanks! I understand the constraint that different fields have very different efficient algorithms and that the standard library should not necessarily include such special case algorithms.

I still have two questions:

Am I correct in assuming that an â€śideologically pureâ€ť implementation of `det` would look like this:

1. It would have a general purpose method that works on `AbstractMatrix{Any}` by using the `*` and `+` operators provided by the field;
2. It would have special fast implementation for `AbstractMatrix{Float}` that uses LU factorization;
3. It would have special fast implementation for `AbstractMatrix{Bool}` that uses something from a number theory package;
etcetera

My second question is: Is there a â€śblessed by the communityâ€ť package for number theory in Julia?

Most Julia programmers try to implement generic algorithms when possible, then provide specialized methods when applicable.

I would not call this an ideology, itâ€™s just a practical approach.

Note that an implementation using `*` and `+` is tricky because, as discussed in the other thread, overflow becomes an issue. `Int` and friends are not \mathbb{Z} but a subset, so it is (technically) not a field.

Actually, the overflow is exactly what I desire. For instance, working with the finite field GF(n) (i.e. Z_n with some extra structure) requires overflow after the n-th element.

I would not call this overflow, itâ€™s just the property of the operator for that field.

The term overflow is usually used for situations where the â€ścorrectâ€ť value (according to the definitions people usually use for the field that is approximated) is outside what you can represent.

If you define the `+` and `*` operators for your (finite) field, and they donâ€™t overflow, there is no reason a generic `det` should not work. If you package the code for this, just include a method for `LinearAlgebra.det` and you should be fine.

It just so happens that for the usual sets people work with (in particular, floats), overflow and precision loss is a major issue, so linear algebra code includes specialized methods for this (and nearly everything else â€” very few calculations will match the formulas usually presented in numerically-naive textbooks).

Nemo.jl contains what you need. Outside of julia, sagemath has pretty good finite field BLAS; I never got around to benchmarking the BLAS performance of nemo.

The issue is that it is extremely difficult to implement an efficient `det` algorithm that does not use division. With division (e.g. via Gaussian elimination), you can do it in O(nÂł) operations. The naive combinatorial determinant formula is O(n!), which is worse than exponential. There are apparently division-free algorithms that require O(n^3.2) operations and similar, but they seem to be of the theoretical genre that are mainly intended for complexity bounds rather than practical implementation (i.e., their constant factors are probably huge, and their implementation looks like it requires a small computer-algebra system).

In consequence, the `det` function in the `LinearAlgebra` standard library implements only a simple and efficient algorithm that requires division.

1 Like

Fascinating! I was only aware of O(n^4) algorithms. These look implementable, but would of course have horrible cancellation properties for anything non-exact.

1 Like