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

use the LinearAlgebraX package for linear algebra? #281

Open
nsajko opened this issue Nov 23, 2023 · 3 comments · May be fixed by #282
Open

use the LinearAlgebraX package for linear algebra? #281

nsajko opened this issue Nov 23, 2023 · 3 comments · May be fixed by #282

Comments

@nsajko
Copy link
Contributor

nsajko commented Nov 23, 2023

The LinearAlgebra.det method defined in this package is not very good:

  1. It's restricted to Matrix, so LinearAlgebra.det fails for, e.g., LinearAlgebra.Symmetric polynomial matrices
  2. It doesn't check that the matrix is square
  3. It uses recursion with depth that's not statically known
  4. It uses sum without an init argument, causing a dynamic dispatch warning with JET

The LinearAlgebraX package, however, provides two different methods that both work with polynomial matrices (as far as I tried): detx and cofactor_det.

Perhaps LinearAlgebraX could be used to improve the linear algebra functionality in this package.

@nsajko
Copy link
Contributor Author

nsajko commented Nov 24, 2023

Slight correction: only cofactor_det works for polynomial matrices, while detx uses it as a fallback in this case, after a try-catch. So I think that LinearAlgebra.det should be implemented using LinearAlgebraX.cofactor_det here.

@nsajko
Copy link
Contributor Author

nsajko commented Nov 24, 2023

Even for Matrix, which is supported by both det and cofactor_det, the latter is much faster. I also tried improving the current code for det, but it was still much slower than cofactor_det for the matrix I tried:

function det_mapreduce_map(M)
    f = i -> Int8(iseven(i) ? -1 : 1)
    let M = M, f = f
        i -> f(i) * M[i, 1] * LinearAlgebra.det(M[1:end.!=i, 2:end])
    end
end

function det_impl(M::Matrix{<:AbstractPolynomialLike})
    m = size(M)[1]
    if m > 2
        mapreduce(
            det_mapreduce_map(M), +, Base.OneTo(m), init = zero(eltype(M)),
        )
    elseif m == 2
        return M[1, 1] * M[2, 2] - M[2, 1] * M[1, 2]
    else
        return M[1, 1]
    end
end

function LinearAlgebra.det(M::Matrix{<:AbstractPolynomialLike})
    LinearAlgebra.checksquare(M)
    det_impl(M)
end

I guess there are two options now:

  1. Depend on LinearAlgebraX
  2. Copy the (short) cofactor_det implementation into this package

What do you think?

@blegat
Copy link
Member

blegat commented Nov 24, 2023

I'm a bit hesitant to add this dependency but I would accept a PR improving the current implementation

nsajko added a commit to nsajko/MultivariatePolynomials.jl that referenced this issue Nov 25, 2023
Still broken for `LinearAlgebra.Symmetric` polynomial matrices,
producing a `MethodError` because of a missing `oneunit` method. This,
however, seems like a separate matter that would better be addressed
by a separate pull request.

Performance comparison:

```
$ ./julia -t 8
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.11.0-DEV.972 (2023-11-23)
 _/ |\__'_|_|_|\__'_|  |  Commit 9884e447e79 (1 day old master)
|__/                   |

julia> using LinearAlgebra, DynamicPolynomials

julia> @PolyVar a b c d e
(a, b, c, d, e)

julia> const m = diagm(
         -2 => fill(a, 14), -1 => fill(b, 15),
         0  => fill(c, 16),
         2  => fill(e, 14), 1  => fill(d, 15))
16×16 Matrix{Polynomial{DynamicPolynomials.Commutative{DynamicPolynomials.CreationOrder}, Graded{LexOrder}, Int64}}:
 c  d  e  0  0  0  0  0  0  0  0  0  0  0  0  0
 b  c  d  e  0  0  0  0  0  0  0  0  0  0  0  0
 a  b  c  d  e  0  0  0  0  0  0  0  0  0  0  0
 0  a  b  c  d  e  0  0  0  0  0  0  0  0  0  0
 0  0  a  b  c  d  e  0  0  0  0  0  0  0  0  0
 0  0  0  a  b  c  d  e  0  0  0  0  0  0  0  0
 0  0  0  0  a  b  c  d  e  0  0  0  0  0  0  0
 0  0  0  0  0  a  b  c  d  e  0  0  0  0  0  0
 0  0  0  0  0  0  a  b  c  d  e  0  0  0  0  0
 0  0  0  0  0  0  0  a  b  c  d  e  0  0  0  0
 0  0  0  0  0  0  0  0  a  b  c  d  e  0  0  0
 0  0  0  0  0  0  0  0  0  a  b  c  d  e  0  0
 0  0  0  0  0  0  0  0  0  0  a  b  c  d  e  0
 0  0  0  0  0  0  0  0  0  0  0  a  b  c  d  e
 0  0  0  0  0  0  0  0  0  0  0  0  a  b  c  d
 0  0  0  0  0  0  0  0  0  0  0  0  0  a  b  c

julia> @time det(m);
  9.301086 seconds (162.91 M allocations: 5.968 GiB, 9.55% gc time, 1.93% compilation time)

julia> @time det(m);
  9.197843 seconds (162.88 M allocations: 5.967 GiB, 10.97% gc time)

julia> @time det(m);
  9.628294 seconds (162.88 M allocations: 5.967 GiB, 10.74% gc time)
```

The above REPL session is with this commit applied. The same
computation with MultivariatePolynomials v0.5.3 ran for multiple
minutes before I decided to just kill it.

Fixes JuliaAlgebra#281
@nsajko nsajko linked a pull request Nov 25, 2023 that will close this issue
nsajko added a commit to nsajko/MultivariatePolynomials.jl that referenced this issue Nov 25, 2023
Still broken for `LinearAlgebra.Symmetric` polynomial matrices,
producing a `MethodError` because of a missing `oneunit` method. This,
however, seems like a separate matter that would better be addressed
by a separate pull request.

Performance comparison:

```
julia> versioninfo()
Julia Version 1.11.0-DEV.972
Commit 9884e447e79 (2023-11-23 16:16 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × AMD Ryzen 3 5300U with Radeon Graphics
  WORD_SIZE: 64
  LLVM: libLLVM-15.0.7 (ORCJIT, znver2)
  Threads: 11 on 8 virtual cores

julia> using LinearAlgebra, DynamicPolynomials

julia> @PolyVar a b c d e;

julia> f(n) = diagm(
          -2 => fill(a, n - 2), -1 => fill(b, n - 1), 0 => fill(c, n),
           2 => fill(e, n - 2),  1 => fill(d, n - 1),
       )
f (generic function with 1 method)

julia> const m15 = f(15);

julia> const m16 = f(16);

julia> @time det(m15);
  2.849449 seconds (47.98 M allocations: 1.832 GiB, 16.50% gc time, 11.14% compilation time)

julia> @time det(m15);
  2.832748 seconds (47.94 M allocations: 1.830 GiB, 22.05% gc time)

julia> @time det(m16);
  6.476343 seconds (112.61 M allocations: 4.299 GiB, 20.43% gc time)

julia> @time det(m16);
  6.737904 seconds (112.61 M allocations: 4.299 GiB, 22.30% gc time)
```

The above REPL session is with this commit applied. The same
computation with MultivariatePolynomials v0.5.3 ran for multiple
minutes before I decided to just kill it.

Fixes JuliaAlgebra#281
nsajko added a commit to nsajko/MultivariatePolynomials.jl that referenced this issue Nov 25, 2023
Performance improvements, support for more types.

Still broken for `LinearAlgebra.Symmetric` polynomial matrices,
producing a `MethodError` because of a missing `oneunit` method. This,
however, seems like a separate matter that would better be addressed
by a separate pull request.

Performance comparison:

```julia-repl
julia> versioninfo()
Julia Version 1.11.0-DEV.972
Commit 9884e447e79 (2023-11-23 16:16 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × AMD Ryzen 3 5300U with Radeon Graphics
  WORD_SIZE: 64
  LLVM: libLLVM-15.0.7 (ORCJIT, znver2)
  Threads: 1 on 8 virtual cores

julia> using LinearAlgebra, DynamicPolynomials

julia> @PolyVar a b c d e;

julia> f(n) = diagm(
          -2 => fill(a, n - 2), -1 => fill(b, n - 1), 0 => fill(c, n),
           2 => fill(e, n - 2),  1 => fill(d, n - 1),
       )
f (generic function with 1 method)

julia> const m15 = f(15);

julia> const m16 = f(16);

julia> @time det(m15);
  2.834004 seconds (47.98 M allocations: 1.832 GiB, 16.66% gc time, 10.61% compilation time)

julia> @time det(m15);
  2.647003 seconds (47.94 M allocations: 1.830 GiB, 18.63% gc time)

julia> @time det(m16);
  6.441887 seconds (112.61 M allocations: 4.299 GiB, 20.25% gc time)

julia> @time det(m16);
  6.390883 seconds (112.61 M allocations: 4.299 GiB, 19.33% gc time)
```

The above REPL session is with this commit applied. The same
computation with MultivariatePolynomials v0.5.3 ran for multiple
minutes before I decided to just kill it.

Fixes JuliaAlgebra#281
nsajko added a commit to nsajko/MultivariatePolynomials.jl that referenced this issue Nov 26, 2023
Performance improvements, support for more types.

Still broken for `LinearAlgebra.Symmetric` polynomial matrices,
producing a `MethodError` because of a missing `oneunit` method. This,
however, seems like a separate matter that would better be addressed
by a separate pull request.

Performance comparison:

```julia-repl
julia> versioninfo()
Julia Version 1.11.0-DEV.972
Commit 9884e447e79 (2023-11-23 16:16 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × AMD Ryzen 3 5300U with Radeon Graphics
  WORD_SIZE: 64
  LLVM: libLLVM-15.0.7 (ORCJIT, znver2)
  Threads: 11 on 8 virtual cores
Environment:
  JULIA_NUM_PRECOMPILE_TASKS = 3
  JULIA_PKG_PRECOMPILE_AUTO = 0

julia> using LinearAlgebra, DynamicPolynomials

julia> function f(n)
         @PolyVar a b c d e
         diagm(
          -2 => fill(a, n - 2), -1 => fill(b, n - 1), 0 => fill(c, n),
           2 => fill(e, n - 2),  1 => fill(d, n - 1),
         )
       end
f (generic function with 1 method)

julia> const m15 = f(15);

julia> const m16 = f(16);

julia> @time det(m15);
  2.489404 seconds (46.04 M allocations: 2.244 GiB, 18.91% gc time, 13.04% compilation time)

julia> @time det(m15);
  2.231880 seconds (45.94 M allocations: 2.238 GiB, 21.50% gc time)

julia> @time det(m16);
  5.362580 seconds (107.70 M allocations: 5.243 GiB, 23.50% gc time)

julia> @time det(m16);
  5.405048 seconds (107.70 M allocations: 5.243 GiB, 23.65% gc time)
```

The above REPL session is with this commit applied. The same
computation with MultivariatePolynomials v0.5.3 ran for multiple
minutes before I decided to just kill it.

Fixes JuliaAlgebra#281
nsajko added a commit to nsajko/MultivariatePolynomials.jl that referenced this issue Nov 27, 2023
Performance improvements, support for more types.

Still broken for `LinearAlgebra.Symmetric` polynomial matrices,
producing a `MethodError` because of a missing `oneunit` method. This,
however, seems like a separate matter that would better be addressed
by a separate pull request.

Performance comparison:

```julia-repl
julia> versioninfo()
Julia Version 1.11.0-DEV.972
Commit 9884e447e79 (2023-11-23 16:16 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × AMD Ryzen 3 5300U with Radeon Graphics
  WORD_SIZE: 64
  LLVM: libLLVM-15.0.7 (ORCJIT, znver2)
  Threads: 11 on 8 virtual cores

julia> using LinearAlgebra, DynamicPolynomials

julia> function f(n)
         @PolyVar a b c d e
         diagm(
          -2 => fill(a, n - 2), -1 => fill(b, n - 1), 0 => fill(c, n),
           2 => fill(e, n - 2),  1 => fill(d, n - 1),
         )
       end
f (generic function with 1 method)

julia> const m15 = f(15);

julia> const m16 = f(16);

julia> @time det(m15);
  1.945673 seconds (45.22 M allocations: 2.261 GiB, 20.60% gc time, 4.02% compilation time)

julia> @time det(m15);
  1.991062 seconds (45.22 M allocations: 2.261 GiB, 23.74% gc time)

julia> @time det(m16);
  4.596664 seconds (106.67 M allocations: 5.324 GiB, 22.65% gc time)

julia> @time det(m16);
  4.648503 seconds (106.67 M allocations: 5.324 GiB, 22.66% gc time)
```

The above REPL session is with this commit applied, and with all other
recent PRs of mine applied, to MultivariatePolynomials.jl,
DynamicPolynomials.jl, and MutableArithmetics.jl. The same computation
with MultivariatePolynomials v0.5.3 ran for multiple minutes before I
decided to just kill it.

Depends on JuliaAlgebra#285.

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

Successfully merging a pull request may close this issue.

2 participants