[WIP] Fixes for SPD manifold and MMatrix #55
Conversation
There was a problem hiding this comment.
Cool, since. I saw how much faster that might get, I am looking forward to being able to use that. For example one of my optimization problems in on the product manifold P(3)^{112x112x50}, i.e. a very large product manifold of 3x3 SPDs, there it might really be beneficial to use MMatrix (the Matlab version on just one vertical slice is for example https://github.com/kellertuer/MVIRT/blob/master/examples/SPD/CaminoSlice.m)
| @@ -215,7 +215,7 @@ denotes the lower triangular matrix with the diagonal multiplied by $\frac{1}{2} | |||
| function exp!(M::MetricManifold{SymmetricPositiveDefinite{N},LogCholeskyMetric}, y, x, v) where N | |||
| (l,w) = spd_to_cholesky(x,v) | |||
| z = exp!(CholeskySpace{N}(),y,l,w) | |||
| mul!(y,z,z') | |||
| y .= z*z' | |||
There was a problem hiding this comment.
Are you sure this is necessary? I thought the old version was more memory efficient?
There was a problem hiding this comment.
In-place multiplication is mul!!
using LinearAlgebra, BenchmarkTools
A = randn(5,5)
S = similar(A)
fun!(z, y) = z .= y*y'
@btime fun!($S, $A);
@btime mul!($S, $A, $A');The first one (fun!) is not allocation-free.
There was a problem hiding this comment.
Well, mul!(y,z,z') gives wrong answer for MMatrix. I've reported it here: JuliaArrays/StaticArrays.jl#697 . In this code z and y is the same matrix, LinearAlgebra handles it correctly but StaticArrays doesn't. Since one additional allocation won't slow this down too much, having a bit slower but correct implementation seems better until StaticArrays is fixed.
There was a problem hiding this comment.
In mul!, the storage variable must not be aliased with any of the factors, as the docstring clearly states. There can't be an in-place multiplication "while you're multiplying". So your issue over at StaticArrays.jl is not an issue: there is simply no guarantee of what happens when C is aliased with either A or B. What happens if the first argument is different from the second and third for MMatrix, as in mul!(y, z, z')?
There was a problem hiding this comment.
I didn't read that docstring and assumed that the mul! here is correct. It looks like mul! only accidentally works on Array here. The difference between mul! for MMatrix and Matrix only happens when there is aliasing. Here y and z are supposed to be the same matrix, it looks like introducing z is just misleading.
There was a problem hiding this comment.
What? What's the type of z? Where does lmul! come from?
There was a problem hiding this comment.
This is really tedious... Seems like there is no triangular multiplication routines for MMatrix.
There was a problem hiding this comment.
Okay, it seems there is no in-place mul! for triangular [S/M]Matrix, only *, see https://github.com/JuliaArrays/StaticArrays.jl/blob/master/src/triangular.jl. So the line above should be okay for those types. It's suboptimal for Base arrays (and larger static matrices that fall back to BLAS), but maybe these are not much used in this context.
There was a problem hiding this comment.
OK,, we can improve this another time.
test/symmetric_positive_definite.jl
Outdated
| exp_log_atol_multiplier = 8.0 | ||
| if T <: MMatrix{3,3,Float64} | ||
| # eigendecomposition of 3x3 SPD matrices from StaticArrays is not very accurate | ||
| exp_log_atol_multiplier = 5.0e7 |
There was a problem hiding this comment.
e7? sounds like really inaccurate.
There was a problem hiding this comment.
I guess StaticArray chooses fast over accurate but that's indeed quite extreme. The algorithm is good enough for Float32 but for Float64 one would expect something more accurate.
| @@ -104,10 +104,10 @@ function test_manifold(M::Manifold, pts::AbstractVector; | |||
| retract!(M, new_pt, pts[1], tv1) | |||
| @test is_manifold_point(M, new_pt) | |||
| for x ∈ pts | |||
| @test isapprox(M, zero_tangent_vector(M, x), log(M, x, x); atol = eps(eltype(x)) * exp_log_atol_multiplier) | |||
| @test isapprox(M, zero_tangent_vector(M, x), inverse_retract(M, x, x); atol = eps(eltype(x)) * exp_log_atol_multiplier) | |||
| @test isapprox(M, x, zero_tangent_vector(M, x), log(M, x, x); atol = eps(eltype(x)) * exp_log_atol_multiplier) | |||
There was a problem hiding this comment.
good that you noticed, indeed x should be added here in order to compare semantically tangent vectors and not points on M.
That also looks like a great application for |
Codecov Report
@@ Coverage Diff @@
## master #55 +/- ##
=======================================
Coverage 90.71% 90.71%
=======================================
Files 17 17
Lines 1012 1012
=======================================
Hits 918 918
Misses 94 94
Continue to review full report at Codecov.
|
I've made some changes trying to get SPD manifolds working with
MMatrix. Let's see how this goes in CI.Ref. JuliaArrays/StaticArrays.jl#696 and JuliaArrays/StaticArrays.jl#697