Charpoly computation improved #804
Conversation
|
Hello and thanks for your contribution :) Would you mind submitting a version of the PR without any automatic formatting changes applied? |
|
Yeah for sure, I apologize for that formatting stuff. |
Codecov Report
@@ Coverage Diff @@
## master #804 +/- ##
=======================================
Coverage 97.14% 97.14%
=======================================
Files 4 4
Lines 315 315
=======================================
Hits 306 306
Misses 9 9 Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here. |
|
No worries, thanks for updating! It seems like github has been thrown off a bit by some of your changes, it shows the entire file as affected by the diff. I think the PR needs a squash, are you familiar with how to perform one? |
No, I'm not familiar with PR squashing, but I'm looking forward to learn about it. |
|
hehe, it's quite simple, but it can be come messy at times. It boils down to combining several different commits so that they appear as a single commit. It can be done from the command-line interface, but I usually use a tool for the purpose, in my case, I use "git lens" in vscode. |
|
Thank you for the explanation, I'm going to squash the PR then. Hopefully it won't mess up. |
stepanoslejsek
force-pushed
the
master
branch
from
d0e0e58
to
a1fd7c8
Compare
|
Hmm, your squash worked (nice ;) but there's still a whole-file diff. Base on this article |
|
Yeah, running this command
I'm not a windows person so I'm not sure how to fix that :/ I suggest googling how to turn that behavior off in your editor |
|
Oh I see now, I'm using linux, let me fix those invisible characters and squash it once again. |
stepanoslejsek
force-pushed
the
master
branch
from
3771180
to
28fd405
Compare
|
It's still there, not exactly sure why, I'll try something else to deal with it. |
|
I think that your editor might be trying too hard to force its settings upon the file by auto formatting. I saw that you had one commit that did away with the line-ending edits, but it had instead inserted tabs instead of spaces for indentation 😃 Generally, these features are best left turned off when working on open-source packages, unless the package happens to follow exactly the same convention that your editor happens to be configured with. Please don't let this discourage you, I do appreciate your PR and I'm sure you'll find the proper workflow soon! |
|
@stepanoslejsek What editor are you using? I found that Emacs + julia-mode + magit usually work, although I have never tried on Linux. |
|
I am using neovim, but I found the solution. It'll be without those invisible characters. I checked it with cat -A command. |
stepanoslejsek
force-pushed
the
master
branch
from
208c759
to
dd1f2d5
Compare
There was a problem hiding this comment.
Overall this looks good, thanks for this PR!
What I'm missing is some form of evaluation whether or not this implementation is actually faster and more accurate than the existing one.
Benchmarking the speed should be straightforward using BenchmarkTools.jl, and benchmarking the accuracy could potentially be done by comparing to the result computed using BigFloat
| @@ -298,18 +299,40 @@ function siso_ss_to_zpk(sys, i, j) | |||
| return z, eigvals(A), k | |||
| end | |||
|
|
|||
| function charpoly(A::AbstractMatrix{T}) where {T} | |||
| N = size(A, 1) | |||
| poly_factors::Vector{T} = vec(ones(T, N+1, 1)) | |||
There was a problem hiding this comment.
| poly_factors::Vector{T} = vec(ones(T, N+1, 1)) | |
| poly_factors = vec(ones(T, N+1)) |
the type assertion is only beneficial if type inference fails
| function charpoly(A::AbstractMatrix{<:Real}) | ||
| Λ = eigvalsnosort(A) | ||
| return prod(roots2real_poly_factors(Λ)) | ||
| t::Matrix{T} = zeros(N, N) # Preallocation |
There was a problem hiding this comment.
| t::Matrix{T} = zeros(N, N) # Preallocation | |
| t = zeros(T, N, N) # Preallocation |
the element type of a zero matrix is specified like this
| Λ = eigvalsnosort(A) | ||
| return prod(roots2real_poly_factors(Λ)) | ||
| t::Matrix{T} = zeros(N, N) # Preallocation | ||
| P::Matrix{T}, H::Matrix{T} = hessenberg(A) |
There was a problem hiding this comment.
| P::Matrix{T}, H::Matrix{T} = hessenberg(A) | |
| P, H = hessenberg(A) |
| @simd for j=N:-1:1 | ||
| @simd for i=1:j | ||
| @simd for k=N-j:-1:1 | ||
| t[k+1, i] = H[i, j]*t[k, j+1]-H[j+1, j]*t[k, i] | ||
| end | ||
| t[1, i] = H[i, j] | ||
| end | ||
| @simd for k=1:N-j |
There was a problem hiding this comment.
| @simd for j=N:-1:1 | |
| @simd for i=1:j | |
| @simd for k=N-j:-1:1 | |
| t[k+1, i] = H[i, j]*t[k, j+1]-H[j+1, j]*t[k, i] | |
| end | |
| t[1, i] = H[i, j] | |
| end | |
| @simd for k=1:N-j | |
| for j=N:-1:1 | |
| for i=1:j | |
| for k=N-j:-1:1 | |
| t[k+1, i] = H[i, j]*t[k, j+1]-H[j+1, j]*t[k, i] | |
| end | |
| t[1, i] = H[i, j] | |
| end | |
| for k=1:N-j |
@simd is unlikely to do anything here since you are not using @inbounds annotations
| end | ||
| end | ||
| i = N-1:-1:0 | ||
| poly_factors[2:N+1, 1] = (-1).^(N.-i).*t[N.-i, 1] |
There was a problem hiding this comment.
| poly_factors[2:N+1, 1] = (-1).^(N.-i).*t[N.-i, 1] | |
| @. poly_factors[2:N+1, 1] = (-1)^(N-i)*t[N-i, 1] |
| @@ -298,18 +299,40 @@ function siso_ss_to_zpk(sys, i, j) | |||
| return z, eigvals(A), k | |||
| end | |||
|
|
|||
| function charpoly(A::AbstractMatrix{T}) where {T} | |||
There was a problem hiding this comment.
| function charpoly(A::AbstractMatrix{T}) where {T} | |
| """ | |
| Implements: "An accurate and efficient algorithm for the computation of the characteristic polynomial of a general square matrix." | |
| """ | |
| function charpoly(A::AbstractMatrix{T}) where {T} |
|
I did a benchmarking on BigFloat matrix using the old and new method and here is the result of benchmarking using @Btime |
Returning reverse array deleted Formatting fixed, handling case when matrix has 0x0 dimensions Commited conversion.jl Invisible characters deleted Truly deleted invisible characters Removed @simd macros, added docstring
stepanoslejsek
force-pushed
the
master
branch
from
cce0213
to
21e35b2
Compare
|
Wonderful, that's a solid 10x improvement :) I think you might have misunderstood slightly what I meant with using BigFloat though, BigFloats are for comparing accuracy, while performance benchmarks are likely better run using the standard |
|
I see, sorry for misunderstanding, I took a random 10x10 Float64 matrix and use the old and new algorithm and pass trough a BigFloat representation and Float64 representation and then took a L2 norm on each coefficient and here is the result.
|
|
That looks amazing :) Thanks for your nice contribution! |
|
May I ask how you performed the |
I ran into same problem with only |
|
Thanks for the hint, OrdinaryDiffEq loads GenericSchur.jl which contains an implementation of |
I took a task for improving computation of characteristic polynomial as it was mentioned in the #TODO note using the provided article. Testing withing runtest.jl was successful.