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

Add support for calculating eigenvalues and eigenvectors #170

Merged
merged 11 commits into from
Oct 26, 2021
1 change: 1 addition & 0 deletions src/Rotations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@ include("principal_value.jl")
include("rodrigues_params.jl")
include("error_maps.jl")
include("rotation_error.jl")
include("eigen.jl")
include("deprecated.jl")

export
Expand Down
42 changes: 42 additions & 0 deletions src/eigen.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,42 @@
# 3D
function LinearAlgebra.eigvals(R::Rotation{3})
θ = rotation_angle(R)
return SVector(exp(-θ*im), exp(θ*im), 1)
end

function LinearAlgebra.eigvecs(R::Rotation{3})
n3 = normalize(rotation_axis(R))
n1 = normalize(perpendicular_vector(n3))
n2 = normalize(n3×n1)
v1 = normalize(n1 + n2*im)
v2 = normalize(n1*im + n2)
v3 = n3
return hcat(v1,v2,v3)
end

function LinearAlgebra.eigen(R::Rotation{3})
return eigen(AngleAxis(R))
end

function LinearAlgebra.eigen(R::AngleAxis)
λs = eigvals(R)
vs = eigvecs(R)
return Eigen(λs, vs)
end


# 2D
function LinearAlgebra.eigvals(R::Rotation{2})
θ = rotation_angle(R)
return SVector(exp(-θ*im), exp(θ*im))
end

function LinearAlgebra.eigvecs(R::Rotation{2})
return @SMatrix [1/√2 1/√2;im/√2 -im/√2]
end

function LinearAlgebra.eigen(R::Rotation{2})
λs = eigvals(R)
vs = eigvecs(R)
Comment on lines +39 to +40
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How much work is repeated here?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Sorry, I didn't quite understand the question.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I meant, some of the calculation going into the eigenvalues is probably being recalculated in eigvecs. Not the way it is currently written but it usual that it falls out that way.

I doubt it has a massive performance impact here.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah, I've updated code and avoid converting to AngleAxis twice.

Before performance update

julia> R = rand(RotMatrix{3})
3×3 RotMatrix3{Float64} with indices SOneTo(3)×SOneTo(3):
  0.196251   -0.948589   0.248325
  0.979358    0.202123  -0.00188397
 -0.0484051   0.243568   0.968675

julia> @benchmark eigen(R)
BenchmarkTools.Trial: 10000 samples with 434 evaluations.
 Range (min  max):  234.707 ns   2.399 μs  ┊ GC (min  max): 0.00%  89.78%
 Time  (median):     238.424 ns              ┊ GC (median):    0.00%
 Time  (mean ± σ):   246.933 ns ± 73.440 ns  ┊ GC (mean ± σ):  1.26% ±  3.77%

  ▃▆█▆▅▃▁▁▁▁             ▁▁▁        ▁                          ▁
  ████████████▇█▇▇▆▆▇▆▆▅█████▆▆▅▇▇▇█████▇▆▆▅▆▆▇▇▇▇▅▅▅▄▅▄▄▃▄▄▃▄ █
  235 ns        Histogram: log(frequency) by time       315 ns <

 Memory estimate: 208 bytes, allocs estimate: 1.

After performance update

julia> R = rand(RotMatrix{3})
3×3 RotMatrix3{Float64} with indices SOneTo(3)×SOneTo(3):
 -0.81819   -0.492907   -0.295987
 -0.287214  -0.0955672   0.953087
 -0.49807    0.864818   -0.0633776

julia> @benchmark eigen(R)
BenchmarkTools.Trial: 10000 samples with 700 evaluations.
 Range (min  max):  179.454 ns  929.644 ns  ┊ GC (min  max): 0.00%  80.02%
 Time  (median):     183.104 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   187.316 ns ±  40.447 ns  ┊ GC (mean ± σ):  1.17% ±  4.40%

   ▃▇██▇▆▄▃▄▄▃▂▁                                                ▂
  ███████████████▆▆▅▆▄▁▃▁▃▁▁▁▁▃▁▁▁▃▁▁▁▁▁▁▃▁▁▃▃▄▄▅▃▄▆▇▇▇▇▇███▇█▇ █
  179 ns        Histogram: log(frequency) by time        232 ns <

 Memory estimate: 208 bytes, allocs estimate: 1.

return Eigen(λs, vs)
end
57 changes: 57 additions & 0 deletions test/eigen.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,57 @@
@testset "Eigen_3D" begin
all_types = (RotMatrix{3}, AngleAxis, RotationVec,
UnitQuaternion, RodriguesParam, MRP,
RotXYZ, RotYZX, RotZXY, RotXZY, RotYXZ, RotZYX,
RotXYX, RotYZY, RotZXZ, RotXZX, RotYXY, RotZYZ,
RotX, RotY, RotZ,
RotXY, RotYZ, RotZX, RotXZ, RotYX, RotZY)
oneaxis_types = (RotX, RotY, RotZ)

@testset "$(T)" for T in all_types, F in (one, rand)
R = F(T)
λs = eigvals(R)
vs = eigvecs(R)
E = eigen(R)
v1 = vs[:,1]
v2 = vs[:,2]
v3 = vs[:,3]
@test R * vs ≈ transpose(λs) .* vs
@test norm(v1) ≈ 1
@test norm(v2) ≈ 1
@test norm(v3) ≈ 1
@test E.values == λs
@test E.vectors == vs
if !(T in oneaxis_types) && VERSION ≥ v"1.2"
# If the rotation angle is in [0°, 180°], then the eigvals will be equal.
# Note that the randomized RotX (and etc.) have rotation angle in [0°, 360°].
# This needs Julia(≥1.2) to get sorted eigenvalues in a canonical order
# See https://github.com/JuliaLang/julia/pull/21598
@test eigvals(R) ≈ eigvals(collect(R))
end
end
end

@testset "Eigen_2D" begin
all_types = (RotMatrix{2}, Angle2d)

@testset "$(T)" for T in all_types, θ in 0.0:0.1:π
R = T(Angle2d(θ))
λs = eigvals(R)
vs = eigvecs(R)
E = eigen(R)
v1 = vs[:,1]
v2 = vs[:,2]
@test R * vs ≈ transpose(λs) .* vs
@test norm(v1) ≈ 1
@test norm(v2) ≈ 1
@test E.values == λs
@test E.vectors == vs
if VERSION ≥ v"1.2"
# If the rotation angle θ is in [0°, 180°], then the eigvals will be equal.
# Note that the randomized RotX (and etc.) have rotation angle in [0°, 360°].
# This needs Julia(≥1.2) to get sorted eigenvalues in a canonical order
# See https://github.com/JuliaLang/julia/pull/21598
@test eigvals(R) ≈ eigvals(collect(R))
end
end
end
1 change: 1 addition & 0 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -28,6 +28,7 @@ include("rodrigues_params.jl")
include("quatmaps.jl")
include("rotation_error.jl")
include("distribution_tests.jl")
include("eigen.jl")
include("deprecated.jl")

include(joinpath(@__DIR__, "..", "perf", "runbenchmarks.jl"))