From 2002ef584ae50b1534452cb210935e3a08cd8a44 Mon Sep 17 00:00:00 2001 From: Ronny Bergmann Date: Mon, 30 Dec 2024 15:20:13 +0100 Subject: [PATCH] Test the new variants. --- src/plans/nonlinear_least_squares_plan.jl | 35 ++--- src/plans/vectorial_plan.jl | 29 ++-- .../test_nonlinear_least_squares_plan.jl | 133 +++++++++++++++++- 3 files changed, 165 insertions(+), 32 deletions(-) diff --git a/src/plans/nonlinear_least_squares_plan.jl b/src/plans/nonlinear_least_squares_plan.jl index c2c9540973..4cf927a92e 100644 --- a/src/plans/nonlinear_least_squares_plan.jl +++ b/src/plans/nonlinear_least_squares_plan.jl @@ -95,8 +95,11 @@ function get_cost( ) where {E<:AbstractEvaluationType,H<:AbstractManifoldHessianObjective} v = 0.0 for i in 1:length(nlso.objective) - v += get_cost(vector_space(1), nlso.smoothing, get_value(nlso.objective, p, i)^2) + v += get_cost( + vector_space(1), nlso.smoothing, abs(get_value(M, nlso.objective, p, i))^2 + ) end + v *= 0.5 return v end function get_cost( @@ -108,8 +111,8 @@ function get_cost( vector_space=Rn, value_cache=get_value(M, nlso.objective, p), ) where {E<:AbstractEvaluationType,H<:AbstractManifoldHessianObjective} - return sum( - get_cost(vector_space(1), nlso.smoothing, value_cache[i]) for + return 0.5 * sum( + get_cost(vector_space(1), nlso.smoothing, abs(value_cache[i])^2) for i in 1:length(value_cache) ) end @@ -125,8 +128,11 @@ function get_cost( ) where {E<:AbstractEvaluationType} v = 0.0 for i in 1:length(nlso.objective) - v += get_cost(vector_space(1), nlso.smoothing, get_value(nlso.objective, p, i)^2, i) + v += get_value( + vector_space(1), nlso.smoothing, abs(get_value(M, nlso.objective, p, i))^2, i + ) end + v *= 0.5 return v end function get_cost( @@ -138,8 +144,8 @@ function get_cost( vector_space=Rn, value_cache=get_value(M, nlso.objective, p), ) where {E<:AbstractEvaluationType} - return sum( - get_cost(vector_space(1), nlso.smoothing, value_cache[i], i) for + return 0.5 * sum( + get_value(vector_space(1), nlso.smoothing, abs(value_cache[i])^2, i) for i in 1:length(value_cache) ) end @@ -191,6 +197,7 @@ function get_jacobian!( J, nlso::NonlinearLeastSquaresObjective{E,AHVF,<:AbstractVectorGradientFunction}, p; + vector_space=Rn, basis::AbstractBasis=get_basis(nlso.objective.jacobian_type), value_cache=get_value(M, nlso.objective, p), ) where {E,AHVF} @@ -226,7 +233,7 @@ _doc_get_residuals_nlso = """ get_residuals(M::AbstractManifold, nlso::NonlinearLeastSquaresObjective, p) get_residuals!(M::AbstractManifold, V, nlso::NonlinearLeastSquaresObjective, p) -Compute the vector of residuals ``s_i(f_i(p))``, ``i=1,…,n`` given the manifold `M`, +Compute the vector of residuals ``f_i(p)``, ``i=1,…,n`` given the manifold `M`, the [`NonlinearLeastSquaresObjective`](@ref) `nlso` and a current point ``p`` on `M`. # Keyword arguments @@ -258,11 +265,10 @@ function get_residuals!( E,<:AbstractVectorFunction{E,<:ComponentVectorialType},H }, p; - vector_space=Rn, kwargs..., ) where {E<:AbstractEvaluationType,H<:AbstractManifoldHessianObjective} for i in 1:length(nlso.objective) - V[i] = get_cost(vector_space(1), nlso.smoothing, get_value(nlso.objective, p, i)^2) + V[i] = get_value(M, nlso.objective, p, i) end return V end @@ -273,11 +279,10 @@ function get_residuals!( E,<:AbstractVectorFunction{E,<:FunctionVectorialType},H }, p; - vector_space=Rn, value_cache=get_value(M, nlso.objective, p), ) where {E<:AbstractEvaluationType,H<:AbstractManifoldHessianObjective} for i in 1:length(value_cache) - V[i] = get_cost(vector_space(1), nlso.smoothing, value_cache[i]) + V[i] = value_cache[i] end return V end @@ -289,13 +294,10 @@ function get_residuals!( E,<:AbstractVectorFunction{E,<:ComponentVectorialType},<:AbstractVectorFunction }, p; - vector_space=Rn, kwargs..., ) where {E<:AbstractEvaluationType} for i in 1:length(nlso.objective) - V[i] = get_cost( - vector_space(1), nlso.smoothing, get_value(nlso.objective, p, i)^2, i - ) + V[i] = get_value(M, nlso.objective, p, i) end return V end @@ -306,11 +308,10 @@ function get_residuals!( E,<:AbstractVectorFunction{E,<:FunctionVectorialType},<:AbstractVectorFunction }, p; - vector_space=Rn, value_cache=get_value(M, nlso.objective, p), ) where {E<:AbstractEvaluationType} for i in 1:length(value_cache) - V[i] = get_cost(vector_space(1), nlso.smoothing, value_cache[i], i) + V[i] = value_cache[i] end return V end diff --git a/src/plans/vectorial_plan.jl b/src/plans/vectorial_plan.jl index 95689d6abc..7168d96d4f 100644 --- a/src/plans/vectorial_plan.jl +++ b/src/plans/vectorial_plan.jl @@ -21,11 +21,14 @@ for example ``g_i(p) ∈ ℝ^m`` or ``\operatorname{grad} g_i(p) ∈ T_p\mathcal # Fields * `basis` an [`AbstractBasis`](@extref `ManifoldsBase.AbstractBasis`) to indicate the default representation. + +# Constructor + CoordinateVectorialType(basis=DefaultOrthonormalBasis()) """ struct CoordinateVectorialType{B<:AbstractBasis} <: AbstractVectorialType basis::B end - +CoordinateVectorialType() = CoordinateVectorialType(DefaultOrthonormalBasis()) """ get_basis(::AbstractVectorialType) @@ -641,7 +644,7 @@ function get_jacobian( n = vgf.range_dimension d = manifold_dimension(M, p) mP = PowerManifold(M, range, vgf.range_dimension) - gradients = zero_vector(mP, fill(p, pM)) + gradients = zero_vector(mP, fill(p, mP)) vgf.jacobian!!(M, gradients, p) # generate the first row to get an eltype c1 = get_coordinates(M, p, gradients[mP, 1], basis) @@ -711,10 +714,8 @@ function get_jacobian!( n = vgf.range_dimension gradients = vgf.jacobian!!(M, p) mP = PowerManifold(M, range, vgf.range_dimension) - println("range ", range, "mP", mP) for i in 1:n c = @view J[i, :] - println("G ", gradients) get_coordinates!(M, c, p, gradients[mP, i], basis) end return J @@ -748,15 +749,21 @@ end # Part II: mutating vgf – allocating jacobian # (a) We have a single gradient function function get_jacobian!( - M::AbstractManifold, JF, vgf::VGF, p; basis::B=DefaultOrthonormalBasis() + M::AbstractManifold, + JF, + vgf::VGF, + p; + basis::B=DefaultOrthonormalBasis(), + range::AbstractPowerRepresentation=get_range(vgf.jacobian_type), ) where { FT, VGF<:AbstractVectorGradientFunction{<:InplaceEvaluation,FT,<:FunctionVectorialType}, B<:AbstractBasis, } - gradients = zero_vector(mP, fill(p, pM)) + mP = PowerManifold(M, range, vgf.range_dimension) + gradients = zero_vector(mP, fill(p, mP)) vgf.jacobian!!(M, gradients, p) - for i in 1:n + for i in 1:(vgf.range_dimension) JF[i, :] .= get_coordinates(M, p, gradients[mP, i], basis) end return JF @@ -772,7 +779,7 @@ function get_jacobian!( ) where { FT,VGF<:AbstractVectorGradientFunction{<:InplaceEvaluation,FT,<:ComponentVectorialType} } - for i in 1:n + for i in 1:(vgf.range_dimension) vgf.jacobian!![i](M, X, p) JF[i, :] .= get_coordinates(M, p, X, basis) end @@ -847,8 +854,7 @@ function get_gradient( ) n = _vgf_index_to_length(i, vgf.range_dimension) pM = PowerManifold(M, range, n) - P = fill(p, pM) - X = zero_vector(pM, P) + X = zero_vector(pM, fill(p, pM)) return get_gradient!(M, X, vgf, p, i, range) end # (c) Special cases where allocations can be skipped @@ -994,8 +1000,7 @@ function get_gradient!( ) where {FT<:AbstractVectorialType} # a type wise safe way to allocate what usually should yield a n-times-d matrix pM = PowerManifold(M, range, vgf.range_dimension...) - P = fill(p, pM) - Y = zero_vector(pM, P) + Y = zero_vector(pM, fill(p, pM)) JF = reshape( get_coordinates(pM, fill(p, pM), Y, vgf.jacobian_type.basis), power_dimensions(pM)..., diff --git a/test/plans/test_nonlinear_least_squares_plan.jl b/test/plans/test_nonlinear_least_squares_plan.jl index 9fa5bec20e..b6999c87d2 100644 --- a/test/plans/test_nonlinear_least_squares_plan.jl +++ b/test/plans/test_nonlinear_least_squares_plan.jl @@ -1,14 +1,141 @@ using Manifolds, Manopt, Test -@testset "Nonlinear lest squares plane" begin - @testset "Test cost/residual/jacobian cases with smoothing" begin end +@testset "Nonlinear lest squares plan" begin + @testset "Test cost/residual/jacobian cases with smoothing" begin + # a simple nlso objetive on R2 + M = Euclidean(2) + d1 = [1, 0] + d2 = [0, 1] + f1(M, x) = norm(x - d1) + f2(M, x) = norm(x - d2) + f(M, x) = [f1(M, x), f2(M, x)] + # Components + f!(M, V, x) = (V .= [f1(M, x), f2(M, x)]) + j1(M, x) = (x - d1) / norm(x - d1) + j1!(M, X, x) = (X .= (x - d1) / norm(x - d1)) + j2(M, x) = (x - d2) / norm(x - d2) + j2!(M, X, x) = (X .= (x - d2) / norm(x - d2)) + # Function + JF(M, x) = [j1(M, x), j2(M, x)] + JF!(M, JF, x) = (JF .= [j1(M, x), j2(M, x)]) + # Jacobi matrix + J(M, x) = cat(j1(M, x), j2(M, x); dims=2) + J!(M, J, x) = (J .= cat(j1(M, x), j2(M, x); dims=2)) + # Smoothing types + s1 = Manopt.smoothing_factory(:Identity) + s2 = Manopt.smoothing_factory((:Identity, 2)) + + # Test all (new) possible combinations of vectorial cost and Jacobian + # (1) [F]unction (Gradient), [C]omponent (Gradients), [J] Coordinate (Jacobian in Basis) + # (2) [a]llocating [i] inplace + # (3) [s] single smoothing [v] vector smoothing + nlsoFas = NonlinearLeastSquaresObjective( + f, JF, 2; jacobian_type=FunctionVectorialType(), smoothing=s1 + ) + nlsoFav = NonlinearLeastSquaresObjective( + f, JF, 2; jacobian_type=FunctionVectorialType(), smoothing=s2 + ) + nlsoFis = NonlinearLeastSquaresObjective( + f!, + JF!, + 2; + evaluation=InplaceEvaluation(), + jacobian_type=FunctionVectorialType(), + smoothing=s1, + ) + nlsoFiv = NonlinearLeastSquaresObjective( + f!, + JF!, + 2; + evaluation=InplaceEvaluation(), + jacobian_type=FunctionVectorialType(), + smoothing=s2, + ) + + nlsoCas = NonlinearLeastSquaresObjective( + [f1, f2], + [j1, j2], + 2; + function_type=ComponentVectorialType(), + jacobian_type=ComponentVectorialType(), + smoothing=s1, + ) + nlsoCav = NonlinearLeastSquaresObjective( + [f1, f2], + [j1, j2], + 2; + function_type=ComponentVectorialType(), + jacobian_type=ComponentVectorialType(), + smoothing=s2, + ) + nlsoCis = NonlinearLeastSquaresObjective( + [f1, f2], + [j1!, j2!], + 2; + function_type=ComponentVectorialType(), + jacobian_type=ComponentVectorialType(), + evaluation=InplaceEvaluation(), + smoothing=s1, + ) + nlsoCiv = NonlinearLeastSquaresObjective( + [f1, f2], + [j1!, j2!], + 2; + function_type=ComponentVectorialType(), + jacobian_type=ComponentVectorialType(), + evaluation=InplaceEvaluation(), + smoothing=s2, + ) + + nlsoJas = NonlinearLeastSquaresObjective(f, J, 2; smoothing=s1) + nlsoJav = NonlinearLeastSquaresObjective(f, J, 2; smoothing=s2) + nlsoJis = NonlinearLeastSquaresObjective( + f!, J!, 2; evaluation=InplaceEvaluation(), smoothing=s1 + ) + nlsoJiv = NonlinearLeastSquaresObjective( + f!, J!, 2; evaluation=InplaceEvaluation(), smoothing=s2 + ) + + p = [0.5, 0.5] + V = [0.0, 0.0] + Vt = [1 / sqrt(2), 1 / sqrt(2)] + G = zeros(2, 2) + Gt = 1 / sqrt(2) .* [-1.0 1.0; 1.0 -1.0] + for nlso in [ + nlsoFas, + nlsoFav, + nlsoFis, + nlsoFiv, + nlsoCas, + nlsoCav, + nlsoCis, + nlsoCiv, + nlsoJas, + nlsoJav, + nlsoJis, + nlsoJiv, + ] + c = get_cost(M, nlso, p) + @test c ≈ 0.5 + fill!(V, 0.0) + get_residuals!(M, V, nlso, p) + @test V == get_residuals(M, nlso, p) + @test V ≈ Vt + @test 0.5 * sum(abs.(V) .^ 2) ≈ c + fill!(G, 0.0) + get_jacobian!(M, G, nlso, p) + @test G == get_jacobian(M, nlso, p) + @test G == Gt + end + end @testset "Smootthing factory" begin s1 = Manopt.smoothing_factory(:Identity) @test s1 isa ManifoldHessianObjective - + @test Manopt.smoothing_factory(s1) === s1 # Passthrough for mhos s2 = Manopt.smoothing_factory((:Identity, 2)) @test s2 isa VectorHessianFunction @test length(s2) == 2 + @test Manopt.smoothing_factory(s2) === s2 # Passthrough for vhfs s3 = Manopt.smoothing_factory((:Identity, 3.0)) @test s3 isa ManifoldHessianObjective