Inverse of a boolean (or integer) matrix vs factorization

Specific Domains Numerics
, ,
1

My understanding is that factorization is almost always better than explicit inversion for a matrix of real (or complex) numbers (because usually you do not actually need the explicit inverse, rather some other operation in which the explicit inverse is involved and which can be done faster through some matrix factorization).

Is that true, in practice, for integer matrices whose inverses are integer matrices (or more generally, for for matrices over a Z/nZ field, e.g. binary matrices)?

If integer matrix factorization is a thing, is it available in Julia? Some other language? I looked through the Nemo.jl documentation but did not really find anything.

2

LU factorizations can be defined for integer matrices, to give rational matrices, and the same approach would work for finite fields (which would give a result in the field). It won’t be especially high performance for integer matrices, since you will likely need to do the computation in Rational{BigInt} to avoid overflow errors, but for finite fields, I think everything might just work out fine. I’m not aware of any implementations, but rolling your own should be pretty easy, and the performance should be pretty good.

QR factorizations won’t work for this case, since orthogonal matrices will usually have irrational entries.

1 Like
3 
function my_lu(A::AbstractMatrix{T}) where {T}
	T <: Integer && return my_lu!(Rational{BigInt}.(A))
	return my_lu!(copy(A))
end

function my_lu!(A::AbstractMatrix{T}) where {T}
    # Extract values
    m, n = size(A)
    minmn = min(m,n)
    info = 0
    @inbounds begin
        for k = 1:minmn
            if !iszero(A[k,k])
                # Scale first column
                if T <: Rational
					Akk = A[k,k]
					for i = k+1:m
						A[i,k] //= Akk
					end
				else
					Akkinv = inv(A[k,k])
					for i = k+1:m
						A[i,k] *= Akkinv
					end
				end
            elseif info == 0
                info = k
            end
            # Update the rest
            for j = k+1:n
                for i = k+1:m
                    A[i,j] -= A[i,k]*A[k,j]
                end
            end
        end
    end
    return LU{T,typeof(A)}(A, 1:minmn, Int(info))
end

Here is a quick implementation that I think works (not heavily tested though)

2 Likes
4

As cool as this is, it’s four orders of magnitudes slower than floating point lu:

julia> bg = BenchmarkGroup()
0-element BenchmarkTools.BenchmarkGroup:
  tags: []

julia> bg[:lu] = @benchmarkable lu(A)\b setup=(A=rand(Int64, 50, 50); b=rand(Int64, 50))
Benchmark(evals=1, seconds=5.0, samples=10000)

julia> bg[:inv] = @benchmarkable inv(A)*b setup=(A=rand(Int64, 50, 50); b=rand(Int64, 50))
Benchmark(evals=1, seconds=5.0, samples=10000)

julia> bg[:my_lu] = @benchmarkable my_lu(A)\b setup=(A=rand(Int64, 50, 50); b=rand(Int64, 50))
Benchmark(evals=1, seconds=5.0, samples=10000)

julia> res = run(bg)
3-element BenchmarkTools.BenchmarkGroup:
  tags: []
  "lu" => Trial(44.042 μs)
  "my_lu" => Trial(533.832 ms)
  "inv" => Trial(60.556 μs)
5

There are many scenarios where it is worth it to wait for the correct answer, rather than getting the wrong one quickly…

There’s actually no need to write a new LU factorisation routine. According to the documentation of LinearAlgebra.lu, it supports any element type that has +, -, *, and /. (If pivoting is chosen (default) the element type should also support abs and <.)

julia> A = rand(Int64, 50, 50);

julia> my_lu(A) == lu(Rational{BigInt}.(A), Val(false))
true
6

You’re right of course. I don’t know why I thought the point would be speed.

7

The factorization and inverse as written in Base already work over finite fields:

julia> struct Mod{p} <: Number
          val::Int
          function Mod{p}(a::Integer) where {p}
                  new(mod(a,p))
          end
       end

julia> Base.zero(::Type{Mod{p}}) where p = Mod{p}(0)

julia> Base.:+(x::Mod{p}, y::Mod{p}) where p = Mod{p}(x.val+y.val)

julia> Base.:*(x::Mod{p}, y::Mod{p}) where p = Mod{p}(x.val*y.val)

julia> Base.:-(x::Mod{p}) where {p} = Mod{p}(-x.val)

julia> Base.:-(x::Mod{p}, y::Mod{p}) where p =x+(-y)

julia> Base.inv(x::Mod{p}) where p = Mod{p}(invmod(x.val,p))

julia> Base.:/(x::Mod{p}, y::Mod{p}) where p = x * inv(y)

julia> Base.isless(x::Mod{p}, y::Mod{p}) where p=x.val<y.val

julia> Base.abs(x::Mod)=x

julia> Base.conj(x::Mod)=x

julia> m=Mod{19}.([1 2 3;3 2 1;1 0 0])
3×3 Array{Mod{19},2}:
 Mod{19}(1)  Mod{19}(2)  Mod{19}(3)
 Mod{19}(3)  Mod{19}(2)  Mod{19}(1)
 Mod{19}(1)  Mod{19}(0)  Mod{19}(0)

julia> inv(m)*m
3×3 Array{Mod{19},2}:
 Mod{19}(1)  Mod{19}(0)  Mod{19}(0)
 Mod{19}(0)  Mod{19}(1)  Mod{19}(0)
 Mod{19}(0)  Mod{19}(0)  Mod{19}(1)
2 Likes
8

I looked through the Nemo.jl documentation but did not really find anything.

It does, but the functionality is documented in AbstractAlgebra.jl:

For solving with integers/rationals, LU factorization is garbage. Fraction-free LU factorization (FFLU) is much better. Solving using p-adic lifting is the best for large matrices (n > 20, say).

Solving using p-adic lifting:

julia> A = MatrixSpace(ZZ, 50, 50)(rand(Int64, 50, 50));

julia> b = MatrixSpace(ZZ, 50, 1)(rand(Int64, 50, 1));

julia> @btime solve_rational(A, b);
  11.291 ms (3303 allocations: 4.51 MiB)

LU and FFLU factorization:

julia> @btime lu(MatrixSpace(QQ, 50, 50)(A));
  776.146 ms (129839 allocations: 8.97 MiB)

julia> @btime fflu(A);
  33.745 ms (9344 allocations: 1.77 MiB)

Oddly enough Nemo doesn’t seem to expose functions for solving giving the precomputed LU or FFLU factorizations though.

For finite fields, the story is basically the same as for floating-point numbers (except that there is no rounding error); standard LU factorization is the way to go.

julia> A = MatrixSpace(GF(19), 500, 500)(rand(Int64, 500, 500));

julia> b = MatrixSpace(GF(19), 500, 1)(rand(Int64, 500, 1));

julia> @btime solve(A, b);
  16.945 ms (2759 allocations: 6.23 MiB)

julia> @btime lu(A);
  18.449 ms (2761 allocations: 8.12 MiB)

I guess Nemo’s lu is slightly slower than its solve here because of some silly overhead; internally, its solve just computes the LU factorization and then solves the triangular systems here.

4 Likes