From f713f4636899c06da81d75524333520552d31f2a Mon Sep 17 00:00:00 2001 From: Joshua Lampert <51029046+JoshuaLampert@users.noreply.github.com> Date: Mon, 21 Oct 2024 17:56:18 +0200 Subject: [PATCH] Add least squares approximation (#93) * allow using different node set for centers * add example and test * add unit tests --- examples/interpolation/least_squares_2d.jl | 28 +++++++++++ src/discretization.jl | 3 +- src/interpolation.jl | 56 ++++++++++++++++------ src/visualization.jl | 6 +-- test/test_examples_interpolation.jl | 11 +++++ test/test_unit.jl | 28 +++++++++++ test/test_util.jl | 25 ++++++++-- 7 files changed, 133 insertions(+), 24 deletions(-) create mode 100644 examples/interpolation/least_squares_2d.jl diff --git a/examples/interpolation/least_squares_2d.jl b/examples/interpolation/least_squares_2d.jl new file mode 100644 index 00000000..cee18fcf --- /dev/null +++ b/examples/interpolation/least_squares_2d.jl @@ -0,0 +1,28 @@ +using KernelInterpolation +using Plots + +# interpolate Franke function +function f(x) + 0.75 * exp(-0.25 * ((9 * x[1] - 2)^2 + (9 * x[2] - 2)^2)) + + 0.75 * exp(-(9 * x[1] + 1)^2 / 49 - (9 * x[2] + 1) / 10) + + 0.5 * exp(-0.25 * ((9 * x[1] - 7)^2 + (9 * x[2] - 3)^2)) - + 0.2 * exp(-(9 * x[1] - 4)^2 - (9 * x[2] - 7)^2) +end + +n = 1089 +nodeset = random_hypercube(n; dim = 2) +values = f.(nodeset) .+ 0.03 * randn(n) +M = 81 +centers = random_hypercube(M; dim = 2) + +kernel = ThinPlateSplineKernel{dim(nodeset)}() +ls = interpolate(nodeset, centers, values, kernel) +itp = interpolate(nodeset, values, kernel) + +N = 40 +many_nodes = homogeneous_hypercube(N; dim = 2) + +p1 = plot(many_nodes, ls, st = :surface, training_nodes = false) +p2 = plot(many_nodes, itp, st = :surface, training_nodes = false) + +plot(p1, p2, layout = (2, 1)) diff --git a/src/discretization.jl b/src/discretization.jl index 060d5696..fc2d0da6 100644 --- a/src/discretization.jl +++ b/src/discretization.jl @@ -53,7 +53,8 @@ function solve_stationary(spatial_discretization::SpatialDiscretization{Dim, Rea # Do not support additional polynomial basis for now xx = polyvars(Dim) ps = monomials(xx, 0:-1) - return Interpolation(kernel, merge(nodeset_inner, nodeset_boundary), c, system_matrix, + nodeset = merge(nodeset_inner, nodeset_boundary) + return Interpolation(kernel, nodeset, nodeset, c, system_matrix, ps, xx) end diff --git a/src/interpolation.jl b/src/interpolation.jl index fa04c187..333c9bb0 100644 --- a/src/interpolation.jl +++ b/src/interpolation.jl @@ -19,6 +19,7 @@ struct Interpolation{Kernel, Dim, RealT, A, Monomials, PolyVars} <: AbstractInterpolation{Kernel, Dim, RealT} kernel::Kernel nodeset::NodeSet{Dim, RealT} + centers::NodeSet{Dim, RealT} c::Vector{RealT} system_matrix::A ps::Monomials @@ -69,7 +70,7 @@ interpolant. See also [`coefficients`](@ref) and [`polynomial_coefficients`](@ref). """ -kernel_coefficients(itp::Interpolation) = itp.c[eachindex(nodeset(itp))] +kernel_coefficients(itp::Interpolation) = itp.c[eachindex(itp.centers)] """ polynomial_coefficients(itp::Interpolation) @@ -79,7 +80,7 @@ interpolant. See also [`coefficients`](@ref) and [`kernel_coefficients`](@ref). """ -polynomial_coefficients(itp::Interpolation) = itp.c[(length(nodeset(itp)) + 1):end] +polynomial_coefficients(itp::Interpolation) = itp.c[(length(itp.centers) + 1):end] """ polynomial_basis(itp::Interpolation) @@ -121,7 +122,7 @@ or collocation is used. system_matrix(itp::Interpolation) = itp.system_matrix @doc raw""" - interpolate(nodeset, values, kernel = GaussKernel{dim(nodeset)}(), m = order(kernel)) + interpolate(nodeset, centers = nodeset, values, kernel = GaussKernel{dim(nodeset)}(), m = order(kernel)) Interpolate the `values` evaluated at the nodes in the `nodeset` to a function using the kernel `kernel` and polynomials up to a order `m` (i.e. degree - 1), i.e., determine the coefficients ``c_j`` and ``d_k`` in the expansion @@ -135,8 +136,10 @@ maximum degree of `m - 1`. If `m = 0`, no polynomial is added. The additional co \sum_{j = 1}^N c_jp_k(x_j) = 0, \quad k = 1,\ldots, Q = \begin{pmatrix}m - 1 + d\\d\end{pmatrix} ``` are enforced. Returns an [`Interpolation`](@ref) object. +If `centers` is provided, the interpolant is a least squares approximation with the centers used for the basis. """ -function interpolate(nodeset::NodeSet{Dim, RealT}, values::Vector{RealT}, +function interpolate(nodeset::NodeSet{Dim, RealT}, + values::Vector{RealT}, kernel = GaussKernel{Dim}(); m = order(kernel)) where {Dim, RealT} @assert dim(kernel) == Dim @@ -149,24 +152,46 @@ function interpolate(nodeset::NodeSet{Dim, RealT}, values::Vector{RealT}, k_matrix = kernel_matrix(nodeset, kernel) p_matrix = polynomial_matrix(nodeset, ps) system_matrix = [k_matrix p_matrix - transpose(p_matrix) zeros(q, q)] + p_matrix' zeros(q, q)] b = [values; zeros(q)] symmetric_system_matrix = Symmetric(system_matrix) c = symmetric_system_matrix \ b - return Interpolation(kernel, nodeset, c, symmetric_system_matrix, ps, xx) + return Interpolation(kernel, nodeset, nodeset, c, symmetric_system_matrix, ps, xx) +end + +# Least squares approximation +function interpolate(nodeset::NodeSet{Dim, RealT}, centers::NodeSet{Dim, RealT}, + values::Vector{RealT}, + kernel = GaussKernel{Dim}(); + m = order(kernel)) where {Dim, RealT} + @assert dim(kernel) == Dim + n = length(nodeset) + @assert length(values) == n + xx = polyvars(Dim) + ps = monomials(xx, 0:(m - 1)) + q = length(ps) + + k_matrix = kernel_matrix(nodeset, centers, kernel) + p_matrix1 = polynomial_matrix(nodeset, ps) + p_matrix2 = polynomial_matrix(centers, ps) + system_matrix = [k_matrix p_matrix1 + p_matrix2' zeros(q, q)] + b = [values; zeros(q)] + c = system_matrix \ b + return Interpolation(kernel, nodeset, centers, c, system_matrix, ps, xx) end # Evaluate interpolant function (itp::Interpolation)(x) s = 0 kernel = interpolation_kernel(itp) - xs = nodeset(itp) + xis = itp.centers c = kernel_coefficients(itp) d = polynomial_coefficients(itp) ps = polynomial_basis(itp) xx = polyvars(itp) for j in eachindex(c) - s += c[j] * kernel(x, xs[j]) + s += c[j] * kernel(x, xis[j]) end for k in eachindex(d) @@ -184,22 +209,22 @@ end function (diff_op_or_pde::Union{AbstractDifferentialOperator, AbstractStationaryEquation})(itp::Interpolation, x) kernel = interpolation_kernel(itp) - xs = nodeset(itp) + xis = itp.centers c = kernel_coefficients(itp) s = zero(eltype(x)) for j in eachindex(c) - s += c[j] * diff_op_or_pde(kernel, x, xs[j]) + s += c[j] * diff_op_or_pde(kernel, x, xis[j]) end return s end function (g::Gradient)(itp::Interpolation, x) kernel = interpolation_kernel(itp) - xs = nodeset(itp) + xis = itp.centers c = kernel_coefficients(itp) s = zero(x) for j in eachindex(c) - s += c[j] * g(kernel, x, xs[j]) + s += c[j] * g(kernel, x, xis[j]) end return s end @@ -223,8 +248,8 @@ function kernel_inner_product(itp1, itp2) @assert kernel == interpolation_kernel(itp2) c_f = kernel_coefficients(itp1) c_g = kernel_coefficients(itp2) - xs = nodeset(itp1) - xis = nodeset(itp2) + xs = itp1.centers + xis = itp2.centers s = 0 for i in eachindex(c_f) for j in eachindex(c_g) @@ -274,7 +299,8 @@ function (titp::TemporalInterpolation)(t) # Do not support additional polynomial basis for now xx = polyvars(dim(semi)) ps = monomials(xx, 0:-1) - itp = Interpolation(kernel, merge(nodeset_inner, nodeset_boundary), c, + nodeset = merge(nodeset_inner, nodeset_boundary) + itp = Interpolation(kernel, nodeset, nodeset, c, semi.cache.mass_matrix, ps, xx) return itp end diff --git a/src/visualization.jl b/src/visualization.jl index 2146cbd4..412622f1 100644 --- a/src/visualization.jl +++ b/src/visualization.jl @@ -106,13 +106,13 @@ end elseif dim(nodeset) == 2 if training_nodes @series begin - x = values_along_dim(itp.nodeset, 1) - y = values_along_dim(itp.nodeset, 2) + x = values_along_dim(itp.centers, 1) + y = values_along_dim(itp.centers, 2) seriestype := :scatter markershape --> :star markersize --> 10 label --> "training nodes" - x, y, itp.(itp.nodeset) + x, y, itp.(itp.centers) end end @series begin diff --git a/test/test_examples_interpolation.jl b/test/test_examples_interpolation.jl index a165b60a..c41b3ca2 100644 --- a/test/test_examples_interpolation.jl +++ b/test/test_examples_interpolation.jl @@ -118,3 +118,14 @@ end @test_include_example(joinpath(EXAMPLES_DIR, "interpolation_2d_condition.jl"), ns=5:10) end + +@testitem "least_squares_2d.jl" setup=[ + Setup, + AdditionalImports, + InterpolationExamples +] begin + @test_include_example(joinpath(EXAMPLES_DIR, "least_squares_2d.jl"), + l2=1.2759520194191292, linf=0.19486087346749836, + l2_ls=0.5375130503454387, linf_ls=0.06810374254243684, + interpolation_test=false, least_square_test=true) +end diff --git a/test/test_unit.jl b/test/test_unit.jl index fadc599c..9a9cb8a2 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -675,9 +675,37 @@ end @test length(polynomial_basis(itp)) == binomial(order(itp) - 1 + dim(nodes), dim(nodes)) @test system_matrix(itp) isa Symmetric + @test size(system_matrix(itp)) == (7, 7) @test isapprox(itp([0.5, 0.5]), 1.0) @test isapprox(kernel_norm(itp), 0.0) + # Least squares approximation + centers = NodeSet([0.0 0.0 + 1.0 0.0 + 0.0 1.0]) + itp = @test_nowarn interpolate(nodes, centers, ff, kernel) + expected_coefficients = [ + 0.0, + 0.0, + 0.0, + 0.0, + 1.0, + 1.0] + coeffs = coefficients(itp) + @test length(coeffs) == length(expected_coefficients) + for i in eachindex(coeffs) + @test isapprox(coeffs[i], expected_coefficients[i], atol = 1e-15) + end + @test order(itp) == order(kernel) + @test length(kernel_coefficients(itp)) == length(itp.centers) + @test length(polynomial_coefficients(itp)) == order(itp) + 1 + @test length(polynomial_basis(itp)) == + binomial(order(itp) - 1 + dim(nodes), dim(nodes)) + @test system_matrix(itp) isa Matrix + @test size(system_matrix(itp)) == (7, 6) + @test isapprox(itp([0.5, 0.5]), 1.0) + @test isapprox(kernel_norm(itp), 0.0, atol = 1e-15) + # 1D interpolation and evaluation nodes = NodeSet(LinRange(0.0, 1.0, 10)) f(x) = sinpi(x[1]) diff --git a/test/test_util.jl b/test/test_util.jl index 036adab5..e0547695 100644 --- a/test/test_util.jl +++ b/test/test_util.jl @@ -1,21 +1,27 @@ """ test_include_example(example; l2=nothing, linf=nothing, + l2_ls=nothing, linf_ls=nothing, atol=1e-12, rtol=sqrt(eps()), - args...) + kwargs...) -Test by calling `include_example(example; parameters...)`. +Test by calling `include_example(example; kwargs...)`. By default, only the absence of error output is checked. """ macro test_include_example(example, args...) local l2 = get_kwarg(args, :l2, nothing) local linf = get_kwarg(args, :linf, nothing) + local l2_ls = get_kwarg(args, :l2_ls, nothing) + local linf_ls = get_kwarg(args, :linf_ls, nothing) local atol = get_kwarg(args, :atol, 1e-12) local rtol = get_kwarg(args, :rtol, sqrt(eps())) + local interpolation_test = get_kwarg(args, :interpolation_test, true) + local least_square_test = get_kwarg(args, :least_square_test, false) local pde_test = get_kwarg(args, :pde_test, false) local kwargs = Pair{Symbol, Any}[] for arg in args if (arg.head == :(=) && - !(arg.args[1] in (:l2, :linf, :atol, :rtol, :pde_test))) + !(arg.args[1] in (:l2, :linf, :l2_ls, :linf_ls, :atol, :rtol, + :interpolation_test, :least_square_test, :pde_test))) push!(kwargs, Pair(arg.args...)) end end @@ -32,14 +38,23 @@ macro test_include_example(example, args...) # assumes `many_nodes` and `values` are defined in the example values_test = itp.(nodeset) # Check interpolation at interpolation nodes - @test isapprox(norm(values .- values_test, Inf), 0; - atol = $atol, rtol = $rtol) + if $interpolation_test + @test isapprox(norm(values .- values_test, Inf), 0; + atol = $atol, rtol = $rtol) + end many_values = f.(many_nodes) many_values_test = itp.(many_nodes) @test isapprox(norm(many_values .- many_values_test), $l2; atol = $atol, rtol = $rtol) @test isapprox(norm(many_values .- many_values_test, Inf), $linf; atol = $atol, rtol = $rtol) + if $least_square_test + many_values_ls = ls.(many_nodes) + @test isapprox(norm(many_values .- many_values_ls), $l2_ls; + atol = $atol, rtol = $rtol) + @test isapprox(norm(many_values .- many_values_ls, Inf), $linf_ls; + atol = $atol, rtol = $rtol) + end else # PDE test # assumes `many_nodes`, `nodes_inner` and `nodeset_boundary` are defined in the example