-
Notifications
You must be signed in to change notification settings - Fork 29
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
Comments
Slight correction: only |
Even for 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:
What do you think? |
I'm a bit hesitant to add this dependency but I would accept a PR improving the current implementation |
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
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
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
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
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.
The
LinearAlgebra.det
method defined in this package is not very good:Matrix
, soLinearAlgebra.det
fails for, e.g.,LinearAlgebra.Symmetric
polynomial matricessum
without aninit
argument, causing a dynamic dispatch warning with JETThe LinearAlgebraX package, however, provides two different methods that both work with polynomial matrices (as far as I tried):
detx
andcofactor_det
.Perhaps LinearAlgebraX could be used to improve the linear algebra functionality in this package.
The text was updated successfully, but these errors were encountered: